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ABSTRACT 


A  new  Boussinesq-type  model  for  surface  water  wave  propagation  in  coastal 
regions  is  derived.  The  model  is  fully  nonlinear  and  accurate  to  0(p4),  where  p, 
is  the  wave  number  nondimensionalized  by  the  water  depth.  As  an  extension  to 
the  second  order  model  proposed  by  Nwogu  (1993),  a  new  dependent  variable 
is  defined  as  a  weighted  average  between  the  velocity  potential  at  two  distinct 
water  depths  to  force  the  model  to  have  a  (4,4)  Pade  approximation  of  the  exact 
linear  dispersion  relationship.  The  fourth  order  polynomial  approximation  for 
the  velocity  potential  vertical  profile  represents  a  great  improvement  over  existing 
0(fi2)  models,  specially  over  the  intermediate  to  deep  water  range.  Nonlinear 
effects  including  generation  of  super  and  subharmonics,  and  amplitude  dispersion 
are  investigated.  A  finite-difference  numerical  scheme  is  developed  for  the  1- 
dimensional  version  of  the  model  for  the  free  surface  displacement  and  a  velocity- 
type  variable.  Several  solitary  wave  solutions  are  studied  and  compared  with 
other  models,  as  well  as  the  solution  to  the  full  problem.  Computations  of  waves 
propagating  over  submerged  bars  are  compared  with  laboratory  measurements 
and  with  results  of  the  fully  nonlinear  0(/i2)  Boussinesq  model  by  Wei  et  all 
(1995).  All  computations  show  that  the  present  model  represents  a  considerable 
improvement  over  the  0(/i2)  model. 
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Chapter  1 


INTRODUCTION 


One  of  the  most  important  tasks  of  coastal  engineers  is  to  accurately  predict 
waves  in  coastal  regions.  Wind  generated  waves  (or  surface  waves)  can  propagate 
from  deep  ocean  to  the  continental  shelf  with  little  energy  loss.  As  the  waves  enter 
the  continental  shelf,  several  transformations  occur  directly  and  indirectly  related 
to  the  fact  that  the  waves  can  now  “feel”  the  bottom.  Deep  water  waves  are  quite 
well  understood  in  terms  of  both  linear  and  nonlinear  wave  theories.  Due  to  the 
fact  only  the  upper  layer  of  the  ocean  (about  half  of  a  wave  length)  is  affected  by 
the  presence  of  the  waves,  in  deep  water  the  treatment  of  the  vertical  dependence 
of  the  flow  is  simpler  than  in  nearshore  areas,  as  it  does  not  depend  on  the  bottom 
topography.  When  the  waves  approach  the  coast,  the  water  depth  is  no  longer 
large  compared  to  the  wave  length,  and  the  waves  experience  effects  such  as  re¬ 
fraction,  diffraction,  reflection,  shoaling,  nonlinear  interactions,  and  finally  reach 
the  surf-zone,  where  wave  breaking  takes  place.  At  this  point  sudden  reduction  in 
the  momentum  flux  of  the  wave  causes  effects  such  as  set-up  (cross-shore  increase 
in  the  free  surface  mean  level  to  balance  the  momentum  flux  reduction),  long¬ 
shore  current  (which,  combined  with  the  fact  that  wave  breaking  puts  sediment 
in  suspension,  strongly  contributes  to  sediment  transport),  among  others.  The 
basic  aim  of  this  dissertation  is  to  contribute  an  accurate  and  reasonably  efficient 
model  capable  of  modeling  waves  from  deep  to  shallow  water  up  to  the  point  of 
(but  excluding)  wave  breaking. 
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The  earliest  nonlinear  model  in  water  of  finite  depth  was  due  to  Airy  (1845), 
known  as  nonlinear  long  wave  theory  or  nonlinear  shallow  water  theory.  Airy’s 
theory  assumes  dispersive  effects  to  be  negligible  (flow  is  depth  independent).  This 
is  a  good  approximation  only  for  very  long  waves  (water  depth  is  small  compared 
to  the  wave  length)  such  as  tsunamis,  infra-gravity  waves,  and  waves  whose  length 
is  of  the  order  of  the  width  of  the  ocean  basin,  in  which  case  Coriolis  effects  are 
also  important.  For  surface  gravity  waves,  shallow  water  theory  is  not  applicable 
in  most  of  the  depth  range  in  coastal  regions. 

The  difficulty  in  modeling  surface  wave  propagation  is  the  fact  that  in 
intermediate  to  shallow  water,  the  vertical  structure  of  the  flow  changes  quickly, 
and  therefore  models  need  either  to  be  3-dimensional,  which  is  computationally 
costly,  or,  somehow,  to  eliminate  the  vertical  dependence,  but  account  for  its 
effects  in  the  2-dimensional  (horizontal)  model.  The  latter  has  been  the  choice 
adopted  in  almost  the  totality  of  the  dispersive  nonlinear  wave  models  in  existence, 
including  the  present  work.  We  now  briefly  review  some  of  these  models. 

1.1  Existing  Boussinesq-Type  Models 

The  earliest  model  which  is  independent  of  the  vertical  coordinate,  but  in¬ 
cludes  weakly  dispersive  and  nonlinear  effects,  was  derived  by  Boussinesq  (1872). 
The  model,  which  was  derived  for  horizontal  bottom  only,  assumes  irrotationality 
of  the  flow,  parabolic  vertical  dependence  of  the  horizontal  velocity  (or  velocity 
potential)  and  linear  vertical  dependence  of  the  vertical  velocity,  also,  it  has  the 
depth  averaged  horizontal  velocity  and  the  free  surface  elevation  as  the  depen¬ 
dent  variables.  As  an  aside,  in  this  work,  we  refer  to  Boussinesq-type  (or  simply 
Boussinesq)  models  as  those  which  assume  irrotationality  and  that  the  velocity 
potential  (or  horizontal  velocity)  has  a  polynomial  vertical  dependence,  leading 
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to  a  set  of  equations  governing  the  free  surface  elevation  and  a  vertically  inde¬ 
pendent  horizontal  velocity-related  variable  (e.g.  depth-averaged  velocity,  total 
mass  flux,  velocity  potential  at  the  bottom,  etc).  Korteweg  and  deVries  (1895) 
used  the  same  assumptions  as  Boussinesq  (1872)  and,  eliminating  the  velocity 
variable,  derived  a  single  equation  (known  as  KdV  equation)  governing  the  free 
surface  elevation  of  weakly  nonlinear  weakly  dispersive  long  waves  propagating  in 
one  direction.  The  KdV  equation  appears  in  many  branches  of  physics  and  has 
received  much  attention,  mostly  due  to  the  fact  that  it  has  analytical  solutions 
in  terms  of  non-element  ary  functions.  In  the  case  of  gravity  waves,  well  known 
solutions  for  the  KdV  equation  are  the  cnoidal  wave  and  the  solitary  wave. 

To  overcome  the  horizontal  bottom  limitation  in  the  original  Boussinesq 
model,  which  prevents  the  model  from  being  a  very  useful  tool  in  coastal  engi¬ 
neering,  Mei  and  Mehaute  (1966)  and  Peregrine  (1966)  derived  Boussinesq  mod¬ 
els  for  variable  depth.  The  two  models  are  similar  in  the  sense  that  both  use  the 
same  asymptotic  assumptions  (weak  nonlinearity  and  dispersiveness),  but  Mei  and 
Mehaute  (1966)  used  the  velocity  at  the  bottom  as  a  dependent  variable,  whereas 
Peregrine  (1966)  used  the  depth-averaged  velocity.  Due  to  its  wide  popularity  in 
the  coastal  engineering  community,  the  model  by  Peregrine  (1966)  is  often  called 
the  standard  Boussinesq  model. 

Although  standard  Boussinesq  models  can  predict  wave  transformation  in 
coastal  regions  with  relative  accuracy,  its  range  of  validity  is  limited  to  fairly 
shallow  water  (McCowan,  1987),  since  its  linear  dispersion  relationship  is  only  a 
polynomial  approximation  of  the  exact  one.  In  order  to  be  applicable  in  deeper 
water,  many  authors  have  suggested  extending  the  validity  of  Boussinesq  mod¬ 
els.  These  “extended  Boussinesq  models”  have  adjustable  rational  polynomial 
approximations  for  the  dispersion  relationship,  a  major  improvement  over  the 
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approximation  resulting  from  the  standard  Boussinesq  models. 

Madsen  et  al.  (1991)  and  Madsen  and  Sqrensen  (1992)  included  higher 
order  terms  with  adjustable  coefficients  into  the  standard  Boussinesq  equations 
for  constant  and  variable  depth,  respectively.  The  addition  of  these  terms  resulted 
in  a  rational  polynomial  expansion  as  the  linear  dispersion  relationship.  The 
coefficients  were  then  adjusted  to  give  better  linear  shoaling  coefficient  and  a 
more  accurate  approximation  to  the  exact  linear  dispersion  relationship,  namely, 
the  (2,2)  Pade  approximant  (Witting,  1984). 

By  defining  the  dependent  variable  as  the  velocity  at  an  arbitrary  depth, 
Nwogu  (1993)  achieved  a  rational  polynomial  approximation  to  the  exact  lin¬ 
ear  dispersion  relationship  without  the  need  to  add  higher  order  terms  to  the 
equations.  Although  the  arbitrary  location  could  be  chosen  to  give  a  (2,2)  Pade 
approximant  as  the  linear  dispersion  relationship,  Nwogu  (1993)  chose  an  alter¬ 
native  value  which  minimized  the  error  in  the  linear  phase  speed  over  some  depth 
range.  Chen  and  Liu  (1995)  derived  a  model  analogous  to  Nwogu’s  but  used  a 
velocity  potential  at  an  arbitrary  depth  as  the  dependent  variable. 

Wei  et  al.  (1995)  used  Nwogu’s  approach  to  derive  a  Boussinesq  model 
(referred  to  henceforth  as  the  WKGS  model)  without  the  weak  nonlinearity  re¬ 
striction.  Numerical  computations  showed  a  major  improvement  over  the  weakly 
nonlinear  model  of  Nwogu,  and  compared  well  with  solitary  wave  solutions  of 
the  full  potential  problem  obtained  with  the  boundary  elements  method  by  Grilli 
et  al.  (1989). 

Schaffer  and  Madsen  (1995)  combined  the  idea  by  Madsen  and  Sqrensen 
(1992)  of  including  higher  order  terms  in  the  equations  with  the  approach  of 
redefining  the  dependent  variable  as  in  Nwogu  (1993),  to  obtain  a  Boussinesq 
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model  which  has  an  extremely  accurate  linear  dispersion  relationship  given  by 
a  (4,4)  Pade  approximant,  despite  having  only  second  order  vertical  polynomial 
approximation  for  the  horizontal  velocity. 

Boussinesq-type  equations  can  also  be  derived  using  a  Hamiltonian  formu¬ 
lation  (Broer,  1974,  Broer,  1975,  Broer  et  al. ,  1976),  where  the  canonical  variables 
are  the  free  surface  elevation  and  the  velocity  potential  evaluated  at  the  free  sur¬ 
face.  The  equations  derived  by  Broer  do  not  produce  a  positive  definite  linear 
dispersion  relationship  and  predicts  negative  phase  speeds  in  deeper  water.  Im¬ 
provements  over  Broer’s  model  to  extend  the  validity  of  the  equations  for  shorter 
waves  were  made  by  Van  der  Veen  and  Wubs  (1993)  whose  dispersion  relation¬ 
ship  is  a  (2,2)  Pade  approximant.  Mooiman  (1991)  uses  a  rational  polynomial 
approximation  with  arbitrary  coefficients  for  the  Hamiltonian  operator,  and  the 
coefficients  are  determined  by  a  minimization  of  the  error  in  the  linear  dispersion 
relationship.  Due  to  the  difficulty  involved  in  computing  the  kinetic  energy  part  of 
the  total  Hamiltonian,  all  Hamiltonian-based  models  to  this  date  are  only  weakly 
nonlinear. 

1.2  Other  Existing  Models 

Although  they  are  the  most  widely  used,  Boussinesq  models  are  not  the 
only  type  of  models  that  eliminate  the  vertical  dependence  from  the  full  problem, 
while  including  dispersive  and  nonlinear  effects. 

Serre  (1953)  derived  a  fully  nonlinear  model  consisting  of  equations  for 
the  horizontal  velocity  and  the  free  surface  elevation.  In  the  derivation,  the  flow 
is  not  assumed  irrotational,  the  horizontal  velocity  is  assumed  constant  over  the 
depth  and  the  vertical  velocity  variation  is  assumed  linear.  Interestingly,  the  fully 
nonlinear  version  of  the  standard  Boussinesq  model  (with  depth  average  velocity) 
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given  in  Mei  (1989)  has  the  exact  same  form  as  Serre’s  equations,  with  the  depth- 
averaged  velocity  replacing  Serre’s  constant-over-depth  velocity.  Exact  solitary 
wave  solutions  of  Serre’s  equations  were  found  by  Su  and  Gardner  (1969)  and 
Seabra-Santos  et  al.  (1987),  and  the  extension  to  variable  depth  and  2-dimensional 
horizontal  coordinates  was  presented  by  Seabra-Santos  et  al.  (1987),  where  the 
transformation  of  solitary  waves  over  a  shelf  was  investigated. 

Green  et  al.  (1974)  and  Green  and  Naghdi  (1976)  derived  a  family  of  mod¬ 
els  which  became  known  as  Green-Naghdi  (or  GN)  models.  GN  is  a  very  general 
theory  whose  only  assumption  is  that  the  vertical  and  horizontal  velocities  are 
polynomials  of  order  N  and  (N  —  1),  respectively.  There  is  no  assumption  re¬ 
garding  irrotationality,  nonlinearity,  or  even  dispersiveness  (not  directly,  at  least). 
With  N  =  1,  GN  recovers  Serre’s  model  (Kirby,  1997).  By  a  different  method, 
Shields  and  Webster  (1988)  derived  GN  equations  and  investigated  several  so¬ 
lutions  including  periodic  and  solitary  wave  properties  for  N  =  1,  N  =  2,  and 
N  =  3,  and  a  variable  depth  numerical  implementation  of  the  unsteady  model 
with  N  =  2.  The  results  agreed  well  with  laboratory  measurements  by  Hansen 
and  Svendsen  (1979).  Other  applications  of  GN  models  can  be  found  in  Demir- 
bilek  and  Webster  (1992)  and  Webster  and  Wehausen  (1995). 

1.3  Present  Model  and  Dissertation  Outline 

Although  existing  Boussinesq  models  have  relatively  accurate  linear  disper¬ 
sion  relationship  and  no  weak  nonlinearity  limitation,  the  approximation  to  the 
vertical  profile  of  the  horizontal  velocity  (or  velocity  potential)  variable  is  only 
a  second  order  polynomial,  which  leads  to  a  poor  representation  of  the  internal 
kinematics  in  intermediate  to  deep  water.  The  inaccuracies  are  even  more  evident 
in  the  vertical  velocity,  which  has  a  linear  depth  dependence.  To  overcome  this 
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problem,  we  propose  a  Boussinesq  model  which  approximates  the  horizontal  veloc¬ 
ity  (or  velocity  potential)  by  a  fourth  order  polynomial  in  the  vertical  coordinate, 
and  consequently  a  third  order  polynomial  for  the  vertical  velocity  is  achieved. 
In  order  to  obtain  a  (4,4)  Pade  approximant  for  the  exact  linear  dispersion  rela¬ 
tionship,  a  new  dependent  variable  is  introduced  as  the  weighted  average  of  the 
velocity  potential  at  2  different  elevations  in  the  water  column.  Since  the  weight 
and  the  locations  at  which  the  new  variable  is  defined  are  arbitrary,  they  can  be 
chosen  so  that  the  (4,4)  Pade  approximant  is  achieved.  The  resulting  model  is 
a  fully  nonlinear  fourth  order  (in  terms  of  vertical  polynomial  approximation  to 
the  velocity  potential)  model  with  (4,4)  Pade  approximant  as  the  linear  dispersion 
relationship.  The  model  does  not  include  wave  breaking  or  bottom  friction  effects. 
The  outline  of  the  dissertation  is  now  presented. 

In  Chapter  2  we  derive  the  model  starting  from  the  full  potential  flow 
boundary  value  problem  for  an  inviscid,  incompressible  fluid.  A  pair  of  coupled 
equations  approximating  the  conservation  of  mass  and  momentum  is  obtained  for 
the  free  surface  elevation  and  the  average  (as  defined  above)  velocity  potential.  A 
three-equation  model  is  then  derived  for  the  free  surface  elevation  and  the  average 
horizontal  velocity  vector. 

In  Chapter  3,  we  analyze  some  analytical  properties  of  the  model,  includ¬ 
ing  linear  dispersion,  where  we  obtain  the  (4,4)  Pade  approximant  for  the  dis¬ 
persion  relationship,  second  order  Stokes-type  interaction  in  a  random  sea,  where 
we  present  the  model’s  ability  to  predict  super  and  sub-harmonics  transfer  co¬ 
efficients,  and  third  order  interactions  in  a  narrow  banded  sea,  where  we  derive 
the  model’s  amplitude  dispersion  coefficient  and  also  its  Schrodinger  equation 
governing  the  propagation  of  the  envelope  of  wave  groups. 

A  numerical  implementation  of  the  1-dimensional  version  of  the  model 
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derived  in  Chapter  2  is  presented  in  Chapter  4.  Wave  generation  inside  the  domain 
and  wave  absorption  by  sponge  layers  to  simulate  radiation  boundary  conditions 
are  also  discussed  and  several  examples  are  shown. 

In  Chapter  5,  we  apply  the  numerical  implementation  of  the  model  to  the 
propagation  of  solitary  waves  over  both  constant  and  variable  depth.  Several 
solitary  wave  properties  are  discussed  and  compared  with  other  models’  results. 

In  Chapter  6,  we  compare  the  model’s  predictions  with  data  sets  from  two 
distinct  laboratory  experiments  of  regular  waves  propagating  over  a  submerged 
bar.  The  results  are  also  compared  with  WKGS  model. 

The  conclusions  and  recommendations  for  future  work  are  presented  in 
Chapter  7. 
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Chapter  2 


DERIVATION  OF  FOURTH  ORDER  FULLY 
NONLINEAR  MODEL 


In  this  chapter  we  derive  a  fully  nonlinear  Boussinesq-type  model  based 
on  a  4th  order  vertical  polynomial  for  the  velocity  potential.  A  set  of  equations 
for  a  velocity-type  variable  is  then  given.  We  assume  the  fluid  is  inviscid  and 
incompressible,  and  the  flow  is  irrotational,  so  that  a  velocity  potential  <f>  exists 
and  the  velocity  field  can  be  written  as: 

u  =  V3&  (2.1) 


where  the  fluid  velocity  vector  u  =  (u,v,w),  and  <j>  are  functions  of  the  spatial 
Cartesian  coordinates  x,y,z  and  time  f,  and  V3  is  the  3  dimensional  gradient 
operator  V3  =  (d/dx,d/dy,d/dz). 


The  full  boundary  value  problem  for  potential  flow  is  given  in  terms  of 
nondimensional  variables  by 


1 


v  +  4>t  +  -$ 


&,  +  p2vV 

4>z  +  /^2v/i  •  v</> 

1 


W  +  iW 

r 

rjt  +  SV(j)  •  — \<f)z 

y2 


=  0; 
=  0; 
=  0; 

=  0; 


— h  <  z  <  Sr) 

(2.2) 

z  —  —h 

(2.3) 

z  =  Sr/ 

(2.4) 

z  =  Sr) 

(2.5) 
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x  and  y  are  the  horizontal  coordinates  scaled  by  a  representative  wave  number 
k0  =  2tt/Io,  where  L0  is  a  wave  length,  z  is  the  vertical  coordinate  starting  at  the 
still  water  level  and  pointing  upwards  and  h  is  the  water  depth,  both  scaled  by  a 
typical  depth  ho-  rj  is  the  water  surface  displacement  scaled  by  a,  representative 
amplitude  a0.  Two  dimensionless  parameters  are  apparent:  5  =  a0/h0  and  p2  = 
(k0h0)2.  Time  t  is  scaled  by  (fco(g'ho)1^2)-1,  and  <f>,  the  velocity  potential,  is  scaled 
by  6h0(gh0)1/2.  g  is  the  acceleration  due  to  the  gravitational  field,  and  V  is  the 
2-dimensional  horizontal  (x,y)  gradient  operator. 

Integrating  (2.2)  over  the  water  column  and  using  (2.3)  and  (2.5),  we  obtain 
a  mass  conservation  equation 


r)t  +  V  •  M  =  0; 


(2.6) 


We  now  proceed  to  derive  model  equations  for  waves  over  an  arbitrary  bottom 
fi(x,y),  and  assuming  S  =  0(1)  and  0(y2)  «  1.  We  assume  an  Nth  degree 
polynomial  approximation  for  cf>  in  the  2  coordinate: 

N 

(t)=  52(nMx,V,t),  (2.7) 

n=0 

where 

C  =  (h  +  z),  (2.8) 

(f>n  are  functions  of  the  horizontal  spatial  coordinates  and  time,  and  are  to  be 
determined.  By  taking  the  limit  of  (2.7)  as  C  -»■  0,  it  is  clear  that  </>0  is  the 
velocity  potential  at  the  bottom  (  =  0.  Substituting  (2.7)  into  (2.3),  we  obtain 
an  expression  for  cf>i  in  terms  of  </>0: 

4>i  =  -g2GS/h  •  V^o,  (2.9) 
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where  G  =  (1  +  /i2|V/i|2)  .  Since  we  are  seeking  an  asymptotic  approximation 

for  </>  in  terms  of  the  parameter  small  parameter  /i2,  it  would  be  consistent  if  we 
expanded  G  in  a  binomial  expansion  around  fj,2  —  0.  However,  we  choose  not  to  do 
this  in  order  to  maintain  the  positive  definiteness  of  this  quantity  as  the  bottom 
slope  becomes  steep.  Substituting  (2.7)  into  (2.2),  and  equating  coefficients  of  like 
powers  of  (  to  zero,  we  obtain  the  following  recursion  formula: 

(n  +  2)(n  +  l)</>n_|_2  +  fJ-2  [(n  +  2)(n  +  l)|Vfi|20n+2  +  (n  +  l)V2fi0n+i 

+  2 (n  +  1  )Vh  ■  V(f)n+i  +  V20n]  =  0  (2.10) 

We  now  use  (2.9)  and  (2.10)  to  obtain  02,  <^>3,  in  terms  of  0O-  The  series  is 
truncated  at  IV  =  4,  yielding: 

0  =  0o  -  f  (gV/i  •  V0oC  +  ^GV20oC2) 

+  H4  {  ^G2V2hVh  ■  V0O  +  GVh  •  V  (GVfi  •  V0O)  C2 
+  |^G2V2/iV20 0  +  ^GVh  ■  V  (GV20o) 

+  ^GV2  (GV/i  •  V0O)]  C3  +  ^<?V2  (GV20o)  C4}  (2-11) 

Commensurate  with  the  extension  of  the  velocity  potential  to  G(/i4),  we  seek  to 
derive  a  set  of  model  equations  having  a  corresponding  dispersion  relationship  in 
the  form  of  a  (4,4)  Pade  approximant  (Witting,  1984)  representing  the  approxi¬ 
mation 

(2'12) 

For  the  case  of  approximations  retaining  terms  to  G(//2),  the  goal  of  obtaining  the 
corresponding  (2,2)  Pade  approximant  is  achieved  by  redefining  the  velocity  poten¬ 
tial  in  terms  of  the  value  of  the  potential  at  an  elevation  za  =  fi[(l+2a)1/2  — 1];  a  = 
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—2/5  and  using  the  resulting  reference  value  <j>a  =  <f>(za )  as  the  dependent  vari¬ 
able;  see  Nwogu  (1993),  Chen  and  Liu  (1995)  and  Kirby  (1997).  This  procedure 
is  not  adequate  for  a  (4,4)  Pade  approximant,  as  shown  in  Appendix  A.  Instead, 
we  define  a  new  dependent  variable 

4>  =  fi<l>a  +  (1  -  (3)(f>b  (2.13) 

where  4>a  and  <f>b  are  the  velocity  potentials  at  elevations  z  =  za  and  z  =  zb, 
and  (3  is  a  weight  parameter.  Relationships  between  these  parameters  to  give 
the  appropriate  dispersion  relationship  will  be  presented  in  Chapter  3.  <j>  may  be 
written  in  terms  of  <j> o  using  (2.11)  yielding: 

4>  =  (AhGVh-  V0o  +  ^5h2GVVo) 

+  n4  |  Bh2  ^G2V2hVh  ■  V<A)  +  GVh  •  V  {GVh  •  V<fo) 


+  Ch 3  [iG2V2/iV20o  +  \GVh  •  V  (GVVo) 

.6  3  v  y 

+  ^GV2  (GVh  •  V</>0)  +  ^Dh4GV2  (GV2^ „)},  (2.14) 

where 

A  =  h/3(k  +  za)  +  (l-/3)(k  +  zb)]  (2.15) 

h 

B  =  —  [/3(/i  +  za)2  +  (1  —  (3){h  +  2t)2j  (2-16) 

C  s  i[/3(/i+*„)3  +  (l-/3)(/i  +  Z»)3]  (2.17) 

D  =  Y[/J(/i  +  2a)1  +  (l-/?)(k  +  ^)<]  (2.18) 


Inverting  (2.14)  gives  a  formula  for  <j>Q  in  terms  of  <f>  which  is  substituted  into 
(2.11),  leading  to  an  approximation  to  the  full  velocity  potential  in  terms  of  <f>: 

t  =  4>  +  l>2  [(Tift  -  0  F,(i)  +  (Bh2  -  C2)  FiM]  +  /  [{Ah  -  C)  F,JM 
+  (Bh2  -  C2)  FS)  +  (Ch3  -  C3)  FS)  +  (Bh3  -  C4)  FS)} ,  (2.19) 
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where 


FS)  =  GVh-Vj) 

W)  =  \gv24> 

F3{4>)  =  Vh  •  V  (AhVh  ■  W$)  +  ivh  •  V  (. Bh2\724> ) 

Fa{4>)  =  -V2  (AhVh  ■  V$)  +  jV2  (Bh2V2$) 

2  4 

•  _  I  v2W/i  •  V<£  -  V/i  •  V  (v/i  •  V<£) 

JiW  =  —^V2hV2<£  —  ~Vh  •  V  (V2^)  —  ^V2  (V/i  •  V^) 

W)  =  -^V2(V2^).  (2.20) 

By  substituting  (2.19)  into  (2.6),  and  neglecting  terms  of  0(/x6)  and  higher,  we 
obtain  the  approximate  conservation  of  mass  equation: 

rjt  —  —V  •  M,  (2.21) 


where 


M  =  H^4>  +  |  (A  —  1)  Fi(<j>)  +  2  ^  F2(<£)  Vfc 


+ 

+ 


(Afe  -  |)  VF,(<«  +  ( Bh 2  -  ¥-)  VF,(*)  j 
li“H  { [(A  -  1)  FaOrt  +  2  (Bh-  y)  F4(« 


V/i 


+  3(’ct2-^F5(«  +  4^D/.s-^F6W 
+  ( Ah  - 1)  vf3(0)  +  (b/>2  -  !pj  VF4(0) 

+  (ch?  -  Spj  VF5(«  +  (oh'  -  Up)  VFe(«)  , 


(2.22) 
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and  H  =  h  +  5r].  Substituting  (2.19)  into  (2.4)  and  neglecting  terms  of  0(/u6)  and 
higher,  we  obtain  the  approximate  Bernoulli  equation  evaluated  at  z  =  St]: 


rj  +  i .t  +  S[(Ah-H)FSt)+(Bh2-H2)FStj 

+  fi4[(Ah-H)F3(h)+(Bh2-H2)FSt) 

+  (< Ch 3  -  H3)  F5(4>t)  +  (. Dh 4  -  H4)  Fe{<j)tj 

+  5-  ||v^|2  +  2V0  •  V  [fi2  {(Ah  -  H)  Ft(4)  +  {Bh2  -  H 2)  F2(]>)} 

+  v4{{Ah-H)FS)+(Bh2-H2)F4(j>)  +  (Ch3-H3)FS ) 

+  ( Dh 4  -  H4)  Fe( })} 

+  v4  |V  {(Ah  -  H)  Fi($)  +  {Bh2  -  H2)  F2(4>)}\2 

+  n2  [FS)  +  2 F2(4>)}2  +  2/r4  [FS)  +  2 HFS)  +  3 H2F5(4>) 

+  AH3FS)}  [FS)  +  2 HF2($)]  =  0.  (2.23) 

The  pair  of  equations  (2.22)  and  (2.23)  form  a  fully  nonlinear  Boussinesq-type 
model  based  on  a  velocity  potential-type  variable,  </>.  We  now  define  a  velocity- 
type  vector: 

»(*,  V,  t)  =  fi  +  (1  -  0)  WU, .  (2.24) 

The  relationship  between  u  and  <j)  can  be  found  by  inverting  the  gradient  of  (2.19) 
and  substituting  into  (2.24),  and  is  given  by 

V<£  =  u  —  /i2V/i  [(A  —  1)  F2\  +  2  (B  —  A)  hF22] 

-  n4Vh  [(A  -  1)  (F41  +  F43)  +  2  (B  -  A)  h  ( F42  +  F44 ) 

+  3(C  -  B)h2F45  +  4(D  -C)h3F46],  (2.25) 


where 


F22(u)  =  ^GV-u 

F4i(u)  =  -\Vh\2[(A-l)Vh-u  +  {B  -  A)hV -u] 

F42{  u)  =  -^\/2h[(A-l)Vh-u  +  {B  -  A)hV-u] 

F43(u)  =  Vh  •  V  (AhVh  ■  u)  +  ^ Vh  •  V  (tfh2V  •  u) 

F44(u)  =  ^V2  (AhVh  •  u)  +  jV2  (Bh2V  •  u) 

m  i 

■  -  V2 hVh  ■  u  -  Vh  ■  V  (' Vh  ■  u) 

F45(u)  =  -^V2hV-u-^Vh-  V  (V  •  u)  -  i  V2  (V/i  •  u) 

F46(u)  =  -lv2(V-u).  (2.26) 


Now  we  substitute  (2.25)  into  the  expression  (2.22)  for  M,  and  into  the  gradient 
of  the  Bernoulli  equation  (2.23).  The  resulting  set  of  evolution  equations  are  the 
approximate  conservation  laws  using  the  velocity-type  variable  u,  and  is  given  by 
(2.21)  with: 


M  =  H^u  +  fi2  {^Ah  —  —  'j  (2VhF22  +  VF21)  +  i^Bh2  —  VF22 
+  n4  \(Ah  -  y)  (2 VhF42  +  VF41  +  2VhF44  +  VF43) 


+  Bh2 - —  )  (VF42  +  3VhF45  +  VF44) 


+  Chd  - 


(4 VhF46  +  VF45)  +  (Dh4  -  ^  VF46  | ,  (2.27) 


for  mass  conservation,  and 

Ut  =  -  V77  -  (|u|2)  +  r4  (77,  ut)  +  r2  (77,  u) ,  (2.28) 

for  momentum  conservation.  U,  Ti,  and  T2  are  given  by 

U  =  u  +  ix2  [(A  —  1)  h  (2 VhF22  +  VF21)  +  (B  -  1)  h2VF22] 
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(2.29) 


+  yx4  [(A  -  1  )h  (2 V/xF42  +  VF41  +  2VhF44  +  VF43) 

+  (F  -  1)  h 2  (VF42  +  3V/xF45  +  VF44) 

+  (C  -  1)  /x3  (4V/xF46  +  VF45)  +  (F  -  1)  h4VF46 

r4  =  yx2V  [xyF2it  +  (2hr)  +  xy2)  F22tj 

+  yx4V  |^?7  (F4it  +  F4 3()  +  ^2 hrj  +  ?y2^  (F42t  +  F44t) 

+  (3 /x2xy  +  3hrj2  +  xy3)  F45t  +  (4 h3r,  +  6 h2V2  +  4 hrf  +  xy4)  F45t]  (2.30) 

r2  =  -  yx25V  {u  •  [(Ah  -  H)  (VF2i  +  2 V/xF22)  +  (F/x2  -  F2)  VF22] 

+  ^(F21  +  2FF22)2} 

-  yx4£V  {u  •  [(Ah  -  H )  (VF41  +  2V/xF42  +  VF43  +  2V/xF44) 

+  (F/x2  -  H2)  (VF42  +  VF44  +  3 V/xF45) 

+  (( Ch 3  -  F3)  (VF45  +  4V/xF46)  +  (F/x4  -  F4)  VF46] 

+  i  |(A/x  -  F))  (VF21  +  2V/xF22)  +  (F/x2  -  F2)  VF42|2 
+  i[(F21  +  2FF22)(F41  +  2FF42 

+  F43  +  2FF44  +  3F2F45  +  4F3F46)]}  (2.31) 

Once  u  is  known,  the  actual  velocity  field  can  be  computed  from  substituting 
(2.25)  into  the  gradient  of  (2.19): 

(u,  V )  =  u  +  yx2  [{Ah  -  C)  (VF21  +  2 V/xF22)  +  (F/x2  -  C2)  VF22] 

+  yx4  [(A/x  -  ()  ( VF41  +  2V/xF42  +  VF43  +  2V/xF44) 

+  (F/x2  -  C2)  (VF42  +  VF44  +  3 V/xF45) 

+  (C/x3  -  C3)  (VF45  +  4VhF46)  +  (j Dh 4  -  C4)  VF46]  (2.32) 

xx;  =  — yx2  [F2i  +  2(F22]  —  yx4  [F44  +  F43 

+  2£  (F42  +  F44)  +  3£2F45  -f  4£3F4e] .  (2.33) 

The  free  parameters  /3:  za,  and  will  be  considered  in  Chapter  3,  when  we  study 
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the  dispersion  properties  of  the  model.  To  obtain  a  set  of  0(/i4)  weakly  disper¬ 
sive,  weakly  nonlinear  equations  with  0(5)  —  0(/i2),  terms  of  0(52fi2J83fj,2,5jj,4, 
82fi4,53fi4,84fj,4,85iJ,4)  should  be  neglected.  The  assumption  0(5)  =  0(//2)  is  used 
in  the  standard  Boussinesq  model.  To  recover  the  WKGS  model,  one  must  ne¬ 
glect  0(fi4)  terms  while  keeping  all  terms  proportional  to  powers  of  5,  and  setting 
f3  —  1,  and  with  Nwogu’s  ajv  being  related  to  A  and  B  by 

A  =  \Ib 

B  =  2aN  +  l  (2.34) 

and  u  being  replaced  by  WKGS’s  ua. 

Equations  (2.21)  with  (2.27),  and  (2.28)  form  the  2-dimensional  version  of  the 
0(fi4)  fully  nonlinear  model  .  A  one-dimensional  version  in  x  can  be  obtained 
by  assuming  u  =  u,  and  V  =  d/dx.  We,  hereafter,  shall  refer  to  the  0(ji4)  fully 
nonlinear  model  as  FN4  and  the  0(/i4)  weakly  nonlinear  model  as  WN4. 
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Chapter  3 


ANALYTICAL  PROPERTIES 

3.1  Linear  Properties 

In  this  section  we  consider  the  flat  bottom  linearized  version  of  the  model 
governing  r]  and  </>  derived  in  Chapter  2,  and  analyze  several  linear  properties  of 
the  model. 

3.1.1  Dispersion 

Neglecting  all  terms  containing  S  and  assuming  a  flat  bottom  (h  =  ho  =  1) 
in  (2.22)  and  (2.23)  we  have  the  following  linear  equations.  Mass  conservation: 

I ft  +  v2^  +  ^(s-i)v2v2^ 

+  T(B2-f-f  +  4)vW  =  0  (31) 

Bernoulli  equation: 

n  +  *  +  y(£-l)V2<A, 

+  £  (b2  -  B  -  J  +  i)  V2V2^,  =  0  (3.2) 
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To  analyze  the  dispersion  properties  of  these  equations,  we  assume  the  following 
general  solution  to  the  equations: 


rj  =  aei{x~ut) 


4>  =  bei{x~ut] 


(3-3) 


where  u  is  the  angular  frequency  nondimensionalized  by  ko(gh0 )1/2,  a  and  b  are 
amplitudes,  and  i  =  y/—l.  Substituting  (3.3)  into  (3.1)  and  (3.2)  we  obtain  the 
linear  dispersion  relationship  for  the  model: 


LO2  = 


1  - 1 1 

(B-i) 

\V2  +  \ 

{B2~  f-f +  i 

i)  ^ 

1  —  2 

(B  -  l)p?  +  1 1 

(3.4) 


The  expression  (3.4)  resembles  the  remarkably  good  (4,4)  Pade  approximant  to 
the  exact  linear  dispersion  relationship  u>2  =  tanh  fi/ ji  (Witting,  1984).  For  (3.4) 
to  be  the  (4,4)  Pade  approximant,  we  set  B  —  1/9  and  D  =  5/189,  and  solve 
(2.16)  and  (2.18)  for  parameters  /?,  za,  and  Zb-  Since  we  have  3  unknowns  and  2 
equations,  there  are  an  infinite  number  of  solutions  that  give  the  desired  values 
of  B  and  D.  However,  an  arbitrary  choice  of  j3  can  give  imaginary  values  of  za 
or  Zb  or  values  lying  outside  of  the  fluid  domain,  causing  these  parameters  to  lack 
physical  significance.  The  relationship  between  /?,  z„,  and  Zb  for  a  (4,4)  Pade 
approximant  is  as  follows: 


Zn  = 


Zb  = 


If _ W 


1/2 


1 567(1  -/?) 

( . M . \ 

\  567(1  -0)/ 


+ 


8 


1/2' 


5670(1  —  (3) 


1/2 


1/2' 


1/2 


-  1 


(3.5) 


Figure  3.1  shows  a  plot  of  the  real  part  of  za  and  Zb  as  given  by  (3.5).  Notice  that 
values  of  0  between  0.018  and  0.467  will  give  both  za  and  Zb  to  be  real  values  lying 
inside  the  water  column.  The  remaining  free  parameter  can  be  chosen  to  improve 
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p 

Figure  3.1:  Values  of  za(f3)  (solid),  Zb(f3)  (dash-dot)  as  a  function  of  weighting 
factor  /?,  corresponding  to  the  (4,4)  Pade  approximant  dispersion 
relation. 

depth  dependent  properties  such  as  linear  shoaling.  This  has  not  been  done  in 
the  present  work,  where  we  arbitrarily  set  (3  =  0.2,  and  substitute  this  value  into 
(3.5)  to  obtain  za  =  —0.4095  and  Zb  =  —0.7726. 

Figure  3.2  shows  comparison  between  the  standard  Boussinesq  theory  (depth 
averaged  velocity),  Nwogu’s  formulation,  and  present  model  with  (4,4)  Pade  ap¬ 
proximant  dispersion  relationship,  of  the  ratio  of  the  phase  speed  with  Airy’s  exact 
linear  solution.  It  is  clear  that  the  present  model  has  improved  linear  dispersion 
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Figure  3.2:  Ratio  of  phase  speed  with  Airy’s  exact  linear  solution.  Standard 
Boussinesq  (dot),  Nwogu’s  (2,2)  Pade  (dash-dot),  Present  (4,4)  Pade 
(dash). 

properties  over  the  already  accurate  Nwogu’s  model  and  closely  reproduces  the 
exact  solution  through  intermediate  to  deep  water.  Similarly,  the  linear  group 
velocity,  defined  as  Cg  —  du/dk  is  shown  in  Figure  3.3  and  the  improvement  over 
Nwogu’s  model  is  even  more  evident. 

Alternatively  to  the  (4,4)  Pade  approximant,  one  could  compute  values  for 
/3,  za,  and  Zb  to  minimize  the  errors  in  the  linear  phase  speed  and  group  velocity 
over  some  depth  range.  This  was  the  procedure  used  by  Nwogu  to  obtain  his 
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optimized  parameter  a  —  —0.39.  However,  the  authors  found  that  the  (4,4)  Pade 
approximant  is  already  sufficiently  accurate  for  any  practical  purposes,  and  any 
parameter  optimization  over  a  normal  water  depth  range  would  result  in  minor 
improvement. 


3.1.2  Internal  Kinematics 

The  internal  kinematics  of  the  present  model  can  be  obtained  from  (2.19). 
For  the  sake  of  comparison  with  the  exact  linear  solution,  we  assume  flat  bottom 
with  h  =  1  and  the  linear  theory  solution  <j>  =  eI_a,t.  The  flat  bottom  one¬ 
dimensional  version  of  (2.19)  is: 

4>  =  +  04, 

+  l,.2  (B2  -  BC  -  \d  +  C4)  (3-6) 


We  define  a  function  f\{z)  as  the  velocity  potential  given  by  (3.6)  normalized  by 
its  value  at  position  z  =  0,  which  is  used  as  a  common  reference: 

l-4[B-(l  +  2)2]  +  4[B2-iJ(l  +  z)2-f  +  M 


/l(z) 


•r\B 


l]  +  4 


B2  —  B 


D 


+  1 
6  '  6 


(3.7) 


The  vertical  velocity  component  w  can  be  obtained  by  differentiating  (3.6)  with 
respect  to  z.  Similarly  to  /x,  a  vertical  velocity  profile  function  can  be  obtained 
by  defining  /^(z)  =  w(z)/w{ 0): 


Mz) 


H2  [(1  +  *)]  +  4  ; 

-B{1  +  z)  + 

^2  +  lY 

-B  +  l 

(3.8) 


Figure  3.4  shows  comparisons  of  fi(z)  between  the  exact  linear  solution  given  by 
cosh[yu(l  +  0)]/  cosh[/i],  Nwogu’s  model  and  the  present  model,  for  various  values  of 
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relative  water  depth  fi.  Notice  that  for  moderately  shallow  water,  the  two  models 
reproduce  the  exact  solution  quite  well;  as  p  increases,  Nwogu’s  model  starts  to 
deviate  strongly,  developing  a  reverse  flow  pattern  near  the  bottom  at  //  «  3, 
while  the  present  model  remains  very  accurate.  Only  at  quite  deep  water  the 
present  model  starts  to  deviate  considerably  from  the  exact  solution  developing 
an  inflection  point  at  fj,  ~  4.24  and  a  reverse  flow  at  /i  ~  6.07. 

Figure  3.5  shows  results  similar  to  Figure  3.4  for  /2(z).  The  equivalent 
exact  linear  solution  is  n  sinh|/i(l  +  z)}/ cosh[/x].  Notice  that  Nwogu’s  model  has  a 
linear  vertical  profile  for  w,  a  poor  representation  in  intermediate  to  deep  water. 
The  present  model  stays  close  to  the  exact  solution  for  a  wide  range  of  /j,.  A  reverse 
vertical  flow  starts  to  appear  at  ^  fa  6.07,  where  the  horizontal  flow  develops  the 
inflection  point.  Finally,  Figure  3.6  shows  the  ratio  to  the  exact  linear  solution 
f3A  =  ^  tanh(ju)  of  the  ratio  between  vertical  and  horizontal  velocities  <j)zl<t>x  at 
z  =  0,  for  the  present  model,  and  Nwogu’s  model.  The  expression  for  /3 

can  be  obtained  by  dividing  the  numerator  of  (3.8)  by  the  numerator  of  (3.7): 

w(z  =  0)  _ _ A*2  +  [~B  +  3] _ 

«(*  =  0)  “  l-f[B-  1]  +  £[£2-£-f +  |] 

The  present  model  agrees  better  with  the  exact  linear  solution  than  Nwogu’s 
model  for  a  wide  depth  range. 

3.2  Nonlinear  Properties 

In  the  previous  sections  we  have  seen  that  the  proposed  model  has  excel¬ 
lent  linear  dispersion  properties  as  well  as  greatly  improved  representation  of  the 
internal  flow  kinematics.  It  is  useful  to  analyze  some  of  the  nonlinear  properties 
of  the  model  by  using  analytical  tools  such  as  Stokes’  type  asymptotic  expansions 
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and  multiple  scales  expansions,  and,  since  these  types  of  analysis  have  been  ex¬ 
tensively  applied  and  studied  for  the  full  boundary  value  potential  problem,  we 
can  have  an  idea  of  how  well  the  nonlinear  version  of  the  present  model  would 
perform  by  comparing  some  of  its  nonlinear  properties  with  those  of  the  full  prob¬ 
lem,  and  also  with  WKGS  and  Nwogu’s  model,  keeping  in  mind  that  a  numerical 
implementation  of  WKGS  model  has  already  been  tested  with  a  good  degree  of 
success. 

In  the  following  sections  we  investigate  0(8)  nonlinear  interactions  in  a  ran¬ 
dom  sea,  and  0(52)  evolution  of  a  narrow  banded  spectrum  wave  train,  governed 
by  Schrodinger  equation. 

The  constant  depth  versions  of  the  evolution  equations  for  <f>  and  rj  are: 

<1,  +  V-jff  |v0+i(B-iff2)vv^ 

+  £  (ff2  -  iff  -  iffff2  +  iff4)  VV2V2#  |  =  0  (3.10) 

v  +  A  +  j  |v^|2  +  y  (ff  -  ff2)  v2it 
+  (ff2  -  iff  -  ffff2  +  iff4)  V2v2* 

i  (b  -  ff2)  v2^ 

+  T(B2-\D-BH2  +  \Hi 
+  iv  |  (B  -  ff2)2  (v2^)  |2  +  iff2  (v2f)2 

+  i  (bH2  -  iff4)  v2iv2v24>  =  0  (3.11) 


+  v^-v 
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3.2.1  Second  Order  Interactions  in  Random  Sea 


We  will  now  look  at  generation  of  super  and  sub-harmonics  by  second  order 
Stokes-type  interactions.  It  is  well  known  that  in  intermediate  and  deep  water  the 
first  nonlinear  correction  of  a  linear  wave  solution  is  a  set  of  bound  waves  (waves 
that  are  forced  by  the  nonlinear  interactions),  also  called  the  superharmonics  (re¬ 
sulting  from  sum-wave  interactions)  and  the  corresponding  subharmonics  (result¬ 
ing  from  difference- wave  interactions)  (Hasselmann,  1962).  These  bound  waves 
are  proportional  to  products  of  the  amplitudes  of  solutions  to  the  linear  equations. 
The  constants  of  proportionality  (which  are  functions  of  the  local  depth)  will  be 
referred  to  as  transfer  coefficients.  Nwogu  (1993)  has  investigated  the  generation 
of  these  bound  waves  in  his  extended  Boussinesq  model  and  found  qualitatively 
reasonable  agreement  with  Stokes’  theory.  Madsen  and  Sprensen  (1993)  have 
found  similar  results.  Kirby  and  Wei  (1994)  extended  Nwogu’s  model  to  full  non¬ 
linearity  and  found  that  the  retention  of  terms  proportional  to  8ji2  (which  are 
neglected  in  Nwogu’s  model  and  the  standard  Boussinesq  model)  is  essential  to  a 
more  accurate  prediction  of  the  transfer  coefficients  to  the  level  of  accuracy  im¬ 
plied  by  the  order  of  retained  dispersive  terms  in  the  original  model  equations. 
Here,  we  derive  the  transfer  coefficients  for  the  present  model  and  compare  to 
results  from  previous  models. 

We  proceed  to  investigate  nonlinear  properties  of  our  model  by  introducing 
the  perturbation  expansion: 

r)  =  ?7o  +  Srji  +  82t}2 

<j>  —  4>o  +  8(j>\  +  82(f>2  (3.12) 

into  (3.10)  and  (3.11),  and  ordering  the  equations  by  powers  of  8.  At  each  order 
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0(Sn )  we  obtain: 


Vnt  +  Li(f)n  —  Fn 

Tin  +  L2<t>nt  =  Gn  (3.13) 

where  L\  and  L2  are  the  linear  operators: 

Li  =  V2  +  y  (b  -  V2' V2 

+  -£(b*-----+1-)  V2V2V2 
4  V  3  6  30/ 

L2  =  l  +  y(Z?-l)V2 

+  t(b2-b-?  +  5)v2v2 

and  the  forcing  terms  are  given  by: 

F0  =  0 
Go  —  0 

Fl  =  -V-(t/0V(/»o)- y(B-l)V-{»7oV(vVo)} 

-  T(B2-B-f+Dv'b»v(v2v^)} 

G,  =  -i(V*,)2  +  ^{2,„V2*„-(B-l)V^„-V(vVo)  +  (V2*)2} 

-  x  {  (!  -  2B)  fW'vV  +  (b2  -  b  -  f  +  i)  V*  ■  V  (v2 v2,„) 

+  i(B  -  1)2V  (v2&)  •  V  (V2*,)  +  2  (b  -  i)  (V2*,)  (v2v2*,)  } 

f2  =  -v-OdVAO-v-OfcV*) 

-  y  [(B  - 1)  [v  •  {„V  (V2*,)}  +  V  .  {,„v  (v2*)}] 

+  V-{>,2v(vVo)}] 

-  £  [(b2  -  b  - 1  +  i)  [v  ■  {„v  (v2v2*,)}  +  V  •  {,„v  (v2vV,)}] 

-  (b -!)*•{,*(***)}] 


(3.14) 

(3.15) 
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2 

G2  =  -V&  Wo  -  y  {-2  [»?i  Wot  +  ifeWi*] 

+  (5-i)[v0o-v(vVi)+v^-v(vVo)] 

+  2W0W1  -  ?7oV2^0(  -  2t?0V(/)0  •  V  (Wo)  +  2770  (Wo)'} 

-  ^  [-2  (tf  -  1)  [771  V2V20ot  +  *7oV2V2<M] 

+  (tf2  -  B  -  |  +  1)  [v (pi  •  V  (v2Wo)  +  V0O  •  V  (v2Wi)] 

+  (B  -  1)2V  (vVi)  •  V  (vVo) 

+  2  (s  -  I)  [(Wi)  (v2Wo)  +  (Wo)  (v2Wi)] 

-  (B  -  l)772V2V20ot  -  2  (5  -  i)  rjoV 4>o  •  V  (v2Wo) 

-  2(B  -  1)770 V  (Wo)  •  V  (Wo)  +  4  -  |)  (Wo)  (V2' Wo)  } 

We  now  assume  the  following  random  sea  as  the  solution  to  the  0(1)  problem: 

Vo  =  ]W"coW;  <£o  =  X^nSiiW;  (3.16) 

n  n 


where  an  and  bn  are  nondimensional  amplitudes  of  the  functions  rjo  and  0o,  ipn  = 
kn  •  x  —  ujnt,  kn  is  the  77- component  wavenumber  vector  nondimensionalized  by 
ko,  x  is  the  horizontal  coordinates  vector  nondimensionalized  by  k0.  u>n  is  the  77- 
component  angular  frequency  nondimensionalized  by  k0(gho )1^2.  Substitution  of 
(3.16)  into  the  0(1)  set  of  equations  gives  a  set  of  n  relationships  between  un  and 

kn  =  |kn|: 


1-1  (B-  +  \  (B>  -  B  -  f  +  1)  fk* 


(3.17) 


We  also  find  a  relationship  between  an  and  bn  given  by: 


oj„ 


~  k  K  a"’ 


(3.18) 
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where 


Kn  =  kn  (l  -  ^-k2n  (b--)  +  ^-k4n  (b2  — 

|  2  V  3/  4  \  3  6  30 


(3.19) 


Following  the  standard  perturbation  technique,  we  substitute  the  0(1)  solution 
(3.16)  into  the  right-hand-side  of  the  0(5)  equations  (3.13)  to  find  the  forcing  of 
the  0(5)  problem.  The  forcings  F  and  G  in  the  mass  and  dynamic  equations 
respectively  are: . 

F  -  7  X/  amal  {^mi  sin(</>/  +  4>m)  +  3~mi  sin(<^  —  <^m)} 


where 


umkf  ±  unk2m  +  (u>i  ±  u;m)(k;  •  km) 


U>lUJm 


(3.20) 


=  ■■■■  -w/wm  (k/  •  km)  +  y2  {u2m  k2m k K\ 

+  ufkfkmKm  +  i(J 3  -  1  )unum(k?  +  fc^)(ki  •  km)  ± 

+  ^  {~^  {B  ~  1)  {ulkikmKm  +  umkmkiKi) 

=F  ^  -  i)  umuik2mkf(kl  +  kj) |  ( 


(3.21) 


Equation  (3.20)  is  identical  to  the  full  Stokes’  theory  result,  except  for  the  approx¬ 
imate  dispersion  relationship.  Assuming  0(y2)  <<  1,  and  using  the  dispersion 
relation  and  binomial  expansions  to  eliminate  Ki  and  Kn,  equation  (3.21)  can  be 
rearranged  to: 


g±  _  ~kf  •  km  +  {owjj  +  u4)  ±  +  6 

m‘  UlUm 


(3.22) 
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which  is,  again,  formally  the  same  as  the  full  Stokes’  theory  result  but  with  the 
approximate  dispersion  relationship. 


The  forced  solution  for  rji  can  be  obtained  by  solving  (3.13)  and  is  given 


by: 


m  =22  a™ai  COS(<&  +  <f>m)  +  uml  COS (<f>i  -  4>m)}  (3.23) 

l  m 


z  k 

4  (^)2  - 


JL±  rp± 


(3.24) 


1  -  fil  1 
2  1 

(®-s) 

i  to2 + i 

;)  (kmi)4 

1  _£ 
x  2 

(B  - 1)(*4)2  +  4  ( 

>2-b-!  +  D 

i  (kii)4 

(3.25) 


kmi  =  |k/  ±  km| ,  =  u>/  ±  (3.26) 

W-mi  are  respectively  the  super  and  subharmonics  transfer  coefficients  of  the 
interaction  between  the  (/,  m)  pair  of  waves.  Figures  3.7  and  3.8  show  comparisons 
of  the  ratio  of  % to  Stokes’  solution,  for  Nwogu’s  model,  WKGS  model,  and  the 
present  model.  Note  that  the  poor  representation  of  these  coefficients  at  small 
/i  in  Nwogu’s  model  is  due  to  the  assumption  of  weak  nonlinearity,  as  discussed 
by  Kirby  and  Wei  (1994).  The  present  model  predicts  superharmonic  amplitudes 
very  accurately  over  a  wide  range  of  water  depths.  The  asymptotic  representation 
of  subharmonic  amplitudes  is  also  more  accurate  than  in  previous  models,  but  the 
new  solution  deviates  more  rapidly  from  the  exact  solution  than  do  the  results  of 
the  previous  models.  Using  the  rearranged  form  of  G^u  given  by  (3.22),  a  much 
better  agreement  is  achieved.  Plots  of  the  rearranged  form  of  using  (3.22) 
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Figure  3.7:  Ratio  of  approximate  superharmonic  transfer  coefficients  to  Stokes’ 
solution.  Stokes’  theory  (solid),  Nwogu  (dot),  WKGS  (dash-dot), 
Present  (dash),  Present  rearranged  (thin  dot,  indistinguishable  from 
exact). 

can  also  be  found  in  figures  3.7  and  3.8,  and  the  lines  are  indistinguishable  from 
the  full  Stokes’  theory  results. 

3.2.2  Third  Order  Interactions  in  Narrow  Banded  Sea 

We  now  extend  our  analysis  to  third  order  interactions  by  deriving  a  cubic 
Schrodinger  equation,  which  governs  the  evolution  of  wave  envelope  resulting  from 
the  propagation  of  a  narrow  banded  spectrum  wave  train,  for  the  present  model 
and  comparing  some  of  its  property  with  the  full  boundary  value  problem,  and  also 
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Figure  3.8:  Ratio  of  approximate  subharmonic  transfer  coefficients  to  Stokes’ 
solution.  Stokes’  theory  (solid),  Nwogu  (dot),  WKGS  (dash-dot), 
Present  (dash),  Present  rearranged  (thin  dot,  indistinguishable  from 
exact). 
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with  the  WKGS  fully  nonlinear  second  order  model.  A  more  detailed  derivation 
of  the  equation  for  the  present  model  in  given  in  Appendix  B.  The  derivation  is 
very  similar  to  the  derivation  for  the  full  boundary  value  problem,  which  can  be 
found  in  Mei  (1989)  and  is  done  using  a  standard  WKB  multiple  scales  approach. 
Here  we  outline  the  derivation. 

A  narrow  banded  wave  train  with  carrier  wave  number  ko  and  angular 
frequency  to  is  assumed  to  be  propagating  mainly  in  x  direction.  At  third  order, 
with  the  expansion  (3.12)  only,  the  perturbation  is  singular,  since  the  forcings  are 
resonant  at  that  order.  In  order  to  avoid  this  problem,  the  independent  variables 
need  to  be  “stretched”,  so  the  time  and  space  variables  are  split  into  “fast”  and 
“slow”  contributions: 

t  =  t'  +  St'  +  S2t'  =  t'  +  Tx  +  T2 

x  =  x'  +  8x'  +  82x'  =  x'  +  X\  +  X2 

y  =  8y'  +  52y,  =  Y1+Y2  (3.27) 

Notice  that  y  has  only  slow  scale  contributions.  We  expand  the  dependent  vari¬ 
ables  as: 

y  =  8rji  +  52y  2  +  83r/3 

cj)  =  S(f>i  +  82<f>2  T  83<f>3  (3.28) 

We  then  substitute  (3.27)  and  (3.28)  into  (3.10)  and  (3.11),  and  order  the 
equations  in  a  manner  analogous  to  what  was  done  in  the  previous  section.  We 
assume  the  solution  to  each  order  to  be  of  the  form 

n+1 

yn  =  E  rlnm(X1,X2,YuY2,Tl,T2)eim<x'-“t '> 

m=— (n+1) 

<£„=  E  <l>mn{XuX2,YuY2,T1,T2)<r'lx'-u*)  (3.29) 
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We  then  seek  an  equation  for  the  wave  envelope,  the  leading  order  wave  amplitude, 
in  X\  and  T)  by  relating  the  coefficients  (amplitudes)  in  (3.29)  of  each  order  to 
the  ones  of  the  previous  order.  After  some  algebra,  the  following  equation  is  found 
before  the  substitution  of  4> and  : 


At2  +  CgAx2  —  -ju"AxlXi  —  iCgAy^Yx  +  i<?i  \A\2  A 
I  (  J.  |  1  A  A  (-1-  +  ^l^2)  J,  1  n  /o  on\ 

+  .^|Ao*+  2a,(l  +  C3^  +  C4M4)  AoTM=°  (3'30) 


Where  A  =  2t]0i  is  the  envelope  amplitude,  u>"  =  d2oj/dk2,  and  expressions  for  C\ , 
C3,  and  C4,  are  given  in  Appendix  B.  The  coefficient  a  1  is  also  given  in  Appendix 
B  for  both  the  present  model  and  the  full  potential  problem.  The  final  equation 
for  the  wave  envelope  A  can  be  obtained  after  some  more  algebra  and  is  given  by 

At  -  -  ^CgAviYi  +  i<r62fi2  |A|2  A  +  71  (r) A  =  0  (3.31) 


Where  r  =  8T\,  and  £  =  Xi  —  CgTy.  Without  the  last  term  (see  Appendix  B)  and 
neglecting  Y\  derivatives,  equation  (3.31)  is  the  nonlinear  cubic  Schrodinger  equa¬ 
tion.  The  last  term  can  be  absorbed  by  defining  the  transformation  A'  =  Aelfyi  dr. 
The  coefficient  a  is  the  sum  of  contributions  from  the  wave-wave  interactions  cq, 
and  wave-current  interactions  cr2,  also  given  in  Appendix  B  for  the  present  model 
as  well  as  for  the  full  potential  problem.  If  the  current  component  (terms  involv¬ 
ing  (f>  10)  of  equation  (3.30)  is  neglected,  and  we  assume  uniform  unidirectional 
propagation  (neglect  all  spatial  derivatives),  equation  (3.30)  can  be  integrated  in 
T2,  and  the  solution  for  A  is  given  by: 

A  =  a0e-,(<TiaoT2)  (3.32) 


where  ao  =  |A|.  The  leading  order  solution  of  77  is,  then: 

r/i  =  a0cos(kx  —  uat )  (3.33) 
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where 

ua  =  to  +  (fy)Wo  (3.34) 

The  coefficient  a i,  therefore,  characterizes  the  amplitude  dispersion  occurring  at 
leading  order  due  to  third  order  wave- wave  interactions.  Figure  3.9  shows  compar¬ 
ison  of  the  ratio  <J\  from  the  present  model  and  from  WKGS  model  to  the  Stokes’ 
solution  to  the  full  problem.  Figure  (3.10)  shows  a  similar  comparison  for  <r2.  In 
both  cases,  the  present  model  appears  to  have  a  better  asymptotic  approximation 
to  the  full  problem,  with  excellent  agreement  in  shallower  water  and  acceptable 
agreement  in  intermediate  to  deep  water. 

Another  important  result  that  can  be  derived  from  a  linear  stability  anal¬ 
ysis  of  (3.31)  is  that  concerning  the  stability  of  a  Stokes’  wave  train  to  sideband 
perturbations  (see  Mei,  1989).  It  is  well  known  that,  according  to  such  analysis, 
a  permanent  form  Stokes’  wave  is  a  stable  solution  of  (3.31)  only  in  the  range 
where  to"  and  cr  =  <7i  +  cr2  have  the  same  signs.  A  comparison  of  the  ratio  of  to"  to 
Stokes’  solution  between  the  present  model  and  WKGS  model  is  shown  in  Figure 
3.11.  Notice  that  the  present  model  has  much  better  agreement  of  to"  with  the 
full  model,  compared  to  WKGS  model,  in  any  depth  up  to  deep  water.  Notice 
that  the  change  in  the  sign  of  to"  in  WKGS  model  at  around  fi  =  3.8  will  cause 
a  change  in  the  stability  of  the  model,  as  long  as  the  sign  of  a  remains  correct, 
which  is  the  case  for  both  Boussinesq  models.  For  the  full  boundary  value  prob¬ 
lem,  there  is  a  change  in  the  stability  at  /i  1.36,  when  a  becomes  positive  ( to " 
remains  negative  at  all  depths).  For  water  deeper  than  this  value  Stokes’  waves 
are  unstable  and  any  sideband  perturbation  will  grow  and  modulate  the  solution. 
A  comparison  of  a  (around  the  zero  crossing)  between  the  present  model,  WKGS 
model  and  Stokes’  solution  to  the  full  problem  is  shown  in  Figure  3.12.  In  the 
present  model,  a  changes  sign  in  shallower  water  than  in  the  full  problem,  the 
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Figure  3.9:  Ratio  of  Schrodinger  equation’s  cubic  term  coefficient  to  full  prob¬ 
lem’s  solution.  Wave-wave  interaction  contribution.  Full  boundary 
value  problem  (solid),  WKGS  (dash-dot),  Present  (dash). 

reverse  being  true  for  WKGS  model.  However,  this  is  hardly  a  problem  since  in 
shallow  to  intermediate  water,  the  time  and  space  scale  in  which  these  modula¬ 
tions  take  place  are  much  longer  than  those  associated  with  the  effects  of  depth 
changes.  Also,  it  is  not  clear  whether  or  not  the  discrepancies  at  higher  values  of 
p  are  only  due  to  a  bad  behavior  of  the  perturbation  expansion. 


39 


2  Stokes 


fcr"  cr 


Chapter  4 


NUMERICAL  IMPLEMENTATION 


In  this  chapter  we  present  the  numerical  implementation  of  the  FN4  (and 
WN4)  model  derived  in  Chapter  2.  The  philosophy  behind  this  scheme  follows 
very  closely  the  one  of  WKGS,  but  extended  to  higher  accuracy  for  consistency 
with  the  higher  accuracy  of  the  model  itself.  The  time  integration  is  done  using  a 
high  order  predictor-corrector  scheme  and  the  spatial  derivatives  are  approximated 
with  high  order  finite  differencing.  The  order  of  accuracy  in  all  the  discretized 
terms  in  the  equations  is  such  that  the  truncation  errors,  which  contain  dispersive- 
type  quantities,  are  always  smaller  than  the  highest-order  dispersive  term  in  the 
equations.  This  was  done  to  assure  that  even  when  using  relatively  large  grid 
spacing,  the  dispersion  introduced  by  the  error  due  to  the  discretization  will  not 
overwhelm  the  dispersive  terms  in  the  equations  themselves.  In  WKGS,  this  is 
accomplished  by  making  the  truncation  errors  of  0(jj.4)  when  combined  with  the 
term  being  discretized,  assuming  kAx  =  0{[x).  Since  the  present  model  contains 
0(/j4)  dispersion,  more  accuracy  in  the  approximate  derivatives  is  needed,  so  that 
the  numerical  truncation  leads  to  errors  of  order  higher  than  0(p4).  Among 
the  several  ways  to  implement  the  boundary  conditions,  we  choose  to  use  fully 
reflective  walls,  with  energy  absorbing  sponge  layers  used  near  the  boundary, 
in  order  to  implement  a  radiation  condition.  This  choice  was  made  due  to  the 
simplicity  and  efficiency  of  this  type  of  boundary  condition.  With  this  kind  of 
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formulation,  however,  it  is  necessary  to  include  some  kind  of  wave  generation 
inside  the  domain.  The  wave  generation  is  implemented  by  introducing  a  source 
function  in  the  mass  conservation  equation  acting  on  a  limited  “source  region” 
conveniently  placed  in  the  domain. 

We  rewrite  the  one  dimensional  equations  with  the  source  function  and  the 
sponge  layers  as  follows: 

r)t  =  —Mx  +  fs(x,t),  (4.1) 

Ut  =  —r)x  -  |  (u2)^  +  Ti  (77,  ut)  +  r2  (77,  u )  -  ufd(x),  (4.2) 

where  M ,  U,  Tj,  T2  are  the  one-dimensional  (in  x)  versions  of  the  quantities 
defined  in  (2.27),  (2.29),  (2.30),  and  (2.31).  Later  in  this  chapter  we  will  discuss 
the  introduction  of  the  two  new  terms  in  the  system:  fs(x ,  t),  the  source  function, 
and  ufd(x),  a  dissipation  term  acting  at  the  sponge  layer(s).  At  this  point,  we 
present  the  finite  difference  formulae  and  solution  method  to  the  equations  above. 

4.1  Discretization  and  Solution  Method 

In  this  section  we  present  the  formulas  used  to  approximate  the  partial 
derivatives  and  the  solution  method  to  the  approximate  equations.  We  discretize 
the  spatial  coordinate  x  by:  Xi  =  iAx,  (i  =  0,1,2 ,  ..IV)  and  time  t  by:  t3  = 
jAt,  (j  =  0,1,2,  ..A^).  There  are  three  basic  steps  in  advancing  the  solution 
by  At  in  time:  (i)  the  right-hand-sides  of  (4.1)  and  (4.2)  are  evaluated,  (ii)  the 
equations  are  integrated  in  time  to  solved  for  77  and  U ,  (iii)  u  is  evaluated  from 
U .  Since  we  use  a  predictor-corrector  integration  method  and  we  have  nonlinear 
terms  containing  time  derivatives  on  the  right-hand-side  of  (4.2),  starting  at  the 
corrector  stage  of  step  (i),  steps  (i)  through  (iii)  are  iterated  until  convergence  is 
attained.  We  present  each  of  these  steps  in  the  following  subsections. 
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4.1.1  Evaluation  of  the  Right-Hand-Sides 


As  we  already  stated,  the  finite-difference  approximations  to  the  spatial 
derivatives  in  the  equations  are  done  is  such  a  way  that  the  truncation  error  in 
each  term  of  the  equations  should  lead  to  errors  of  order  higher  than  0(/i4).  Notice 
that  the  momentum  equation  (4.2)  contains  first  order  time  derivatives  of  u  in  I\. 
These  time  derivatives  are  evaluated  in  conjunction  with  the  predictor- corrector 
iterated  scheme  (presented  later)  by  using  a  finite  difference  approximation  with 
values  of  u  at  times  (j,j  —  l,j  —  2 ,j  —  3 ,j  —  4) A t  in  the  predictor  stage  and 
(j  + 1>  3i  j  ~  1>  3  ~ 2,  j  —  3,  j  —  4) At  in  the  corrector  stage.  The  formulas  for  the  time 
derivatives  at  each  of  those  j  locations  were  obtained  by  expanding  the  variables  in 
Taylor  series  around  each  j,  multiplying  each  expansion  by  a  coefficient  and  solving 
the  system  of  equations  resulting  from  setting  the  combination  of  coefficients  of 
the  higher  derivatives  of  t  to  zero  (this  is  the  standard  procedure  to  find  finite 
difference  formulas,  and  was  used  throughout  this  chapter).  The  formulas  are: 

u  ,r4  =  i^(-3fi|  +  16ir‘-36Sr2  +  48ar3-25ir‘)+0(Af‘) 

a,r3  =  — ^  (a?  -  ear1  +  isar2  -  ioar3  -  3ar4)  +  0(a<4) 

a  ,f2  =  (-a;  +  sap1  -  sar3  +  ar ')+  O(AC) 

s,rl  =  (3fi|  +  ioar1  -  isar2  +  ear3  -  iar4)  +  o(At4) 

a,  J  =  (25a;  -  48ar‘  +  sear2  -  iear3  +  3ar4)  +  o(a<4) 

for  the  predictor  stage,  and 

ut4  =  ^  (l2uf+1  -  7bui  +  200uf1  -  300uf2 

+  300ur3  -  137«r4)  +  0(At5) 

fi,r3  =  ^(-3uri+20uJ-606r1  +  120ar2 

-  65af-3  -  i2ar4)  +  0(A(5) 
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Uti  2  =  (2uffl  -  15 u\  +  60 u3i  1  -  20uf  2 

-  30  uf3  +  3uf4)  +  0(Af5) 

«ii_1  =  — ^  (-3u-+1  +  30^  +  20uf1  -  eOuf2 
+  15u/-3  —  2u{~ 4)  +  0(At5) 
uti  =  (12  uj+1  +  65  u\  -  mar1  +  60uf2 

-  20uf3  +  3«r4)  +  0{At5) 

u4+1  =  g^  (l37^+1  -  mui  +  SOOuf1  -  200uf2 

+  75  «r3  -  mr4) + c>(Ai5) 


for  the  corrector  stage. 

The  formulas  for  all  spatial  derivatives  in  (4.1)  and  (4.2)  at  location  iAx 
and  at  time  jAt  are  given  below  for  centered  and  off-centered  (near  but  exclud¬ 
ing  the  boundaries,  where  boundary  conditions  are  used)  locations.  We  use  the 
variable  u  as  an  example. 

For  derivatives  appearing  in  terms  of  0(1)  we  have: 
centered  at  i : 

=  g^  [45  (u{+1  -  uti)  -  9  (u{+2  -  uQ 

+  UJi+ 3  -  Uj_3 ]  +  0( AX6) 

u*xi  =  180^t2  [270  («*+1  +  ar,)  -  27  (u3i+2  +  uj_2) 

+  2  (u-+3  +  tt;_3)  -  490U’]  +  0(Ax5) 


off-centered  at  i  =  1  and  i  —  N  —  1: 


u 


3 

xi,N-l 


± 


1 


60Aa: 


j^— IOuq  ,n  77«{iW_i  +  lSOu^jv— 2 


lOOu 


3 

3,7V— 3 
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+  50u4^_4  —  +  2^6, AT- 6  +  0(Ax 6) 

Uxxi'N-i  =  180Aa,2  [l37 u30N  —  IA7u{'N_x  —  255u32N_2  +  470t^>Ar_3 

—  285u{tN_4  +  93u^n_5  —  13tt^jv_6  +  0(Ax  ) 

off-centered  at  i  =  2  and  i  =  N  ~  2: 

^an,jV  — 1  =  ±gQA;c  ^o,n  ~  —  35£t2,AT— 2  4"  80u3jy_3 

—  30u^iAr_4  +  ~~  ^6,iV-6  4"  0(Ax  ) 

=  180Ax2  — 13ujiW  +  228£tj)Ar_j  —  420u^  ^_2  +  200u^iAr_3 
4"  1^4,A^-4  “  ^2u5  ^_5  +  2 uj6  N_6  +  0(Aa:  ) 

The  “±”  indicates  that  the  expression  should  be  positive  for  the  points  near  the 
left  boundary  ( i  =  0)  and  negative  near  the  right  boundary  (i  -  N).  This 
convention  applies  for  the  formulas  below  as  well. 

For  terms  of  0(/j,2): 

centered  at  i: 

Uxi  =  [8  (Cl  -  «ti)  -  fi?+2  +  <t.]  +  0(Ax4) 

Uxxi  =  12Aj.2  [16  («i+l  +  «i- 1)  “  Si+2  ~  «i-2  “  30^]  +  0(Ax3) 

Uxxxi  =  [-13  (u*+1  -  uf_j)  +  8  (u4+2  -  u\_2) 

-  (^i+3  -  Wi-3)  -  490uf]  +  0(Ax4) 

off-centered  at  i  =  1  and  2  =  iV  —  1: 

=  i  1^2,N-2 

—  &u{n_3  +  w^}Ar_4j  +  0(Ax4) 

Uxxi,N-l  =  [n«o,JV  -  20«1,AT— 1  +  6uJ2,Ar-2 
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+  4utl,N-3  ~  «4,JV-4]  +  0(Ax3) 

Uxxx l,N-l  —  [— 15^0, AT  "I"  56u{tN_  1  —  83u{  n_2  +  64uJ3  N_3 

—  29u34N_4  +  8u(n_3  —  ^6^-6  +  0(Ax4) 

off-centered  at  i  =  2  and  i  =  N  —  2: 

Wx*a:ltJV_j  =  i  8^2-3  ~™0,N  _  ^l,jV-l  d"  ^^2,N-2  ~~  ^^3,N-3 

d*  29 u4  N_4  —  8u3  N_3  4-  ^6,Af~$]  d"  ^(Aa;  ) 

For  terms  of  0(fi4): 

centered  at  i:  , 

Uxi  =  2^  (Si+i  -  “i-i)  +  0(Ax2) 

<W  =  2^2  (“i+1  ~  2ili  +  «<-l)  +  O(Aa0 

=  “3  [-2  ( u3i+1  -  uf_i)  +  u?+2  -  w-_2]  +  0(Ax2) 

=  — ~4  [~4  (“i+1  +«f-i)  d-“i+2  d-«d_2  +6uj]  +0(Aa;) 

off-centered  at  i  =  1  and  i  =  N  —  1: 

=  "^i2Ax^  ~^o,n  d*  10«1at-i  —  2 

d-  6u^_3  -  uj^]  +  0(Az2) 

Wxxxx^AT-X  =  12Aar4  _  d^l,AT-l  d-  6u2iN-2 

-  4^3, AT— 3  +  «4,JV-4]  +  O(Ax) 

For  equation  (4.1)  all  but  the  boundary  points  (that  is,  i  —  1, ...,  IV  —  1)  are 
evaluated.  After  the  time  integration  is  done  (see  next  subsection),  a  boundary 
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condition  T]x  =  0  is  applied  at  i  =  0  and  i  —  N.  A  7-point  off-centered  derivative 
of  rj  is  used  and  t)q+1  and  can  be  obtained  as: 

(360^-,  -  450t$_2  +  400<U_3 

225<J}_4  +  72 -  10^j;_6)  +  0(  Ax7) 

For  equation  (4.2),  we  only  evaluate  the  terms  at  points  i  =  2, ...,  JV  —  2.  The 
remaining  points  do  not  need  to  be  evaluated  since  at  those  points,  the  values  of 
u  are  determined  by  boundary  conditions.  This  is  done  when  we  evaluate  u  from 
U,  defined  in  (2.29),  and  the  procedure  will  be  explained  later. 

4.1.2  Time  Integration 

The  integration  method  used  is  a  5th  order  predictor,  6th  order  corrector, 
Adams-Bashforth-Moulton  scheme.  Once  the  dependent  variables  are  known  at 
times  ( j  —  4,j  —  3 ,j  —  2 ,j  —  l,j)At,  and  the  right-hand-sides  of  the  equations 
have  been  evaluated,  estimates  of  both  rj  and  U  at  time  (j  +  l)Af  are  made  using 
the  predictor  stage: 

v!+lp  =  »?  +  ^  (won?  -  2774V?-1 

+  2616V;.i“2  -  1274 v;j-3  +  25iv;j~4)  +  0(Ate), 

where  index  p  stands  for  predictor,  vj  is  either  U  or  r/,  and  Vf  is  the  right-hand- 
side  of  the  respective  equation,  at  x  =  iAx  and  t  =  jAt.  With  the  estimate  [7J+1p, 
we  evaluate  uJ+1p  (see  next  subsection)  then  estimate  the  right-hand-sides  of  the 
equations  at  V/+1  ,  and  iterate  the  corrector  stage: 

vi  +  Aj  +  1427V;-'  -  798V?-1 

482 V/“2  -  17ZV/~3  +  27V;i-4)  +  0(At7), 


,J+ 1 


+ 
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until  the  error  between  vj+1  and  vf+1p  (where,  again,  v  applies  here  applies  to 
both  rj  and  U )  is  small.  We  define  an  error  estimate  for  the  iteration  process  as: 


Ei 


iter 


N 

E 

i=0 


(1fMp  - 


(’ll 


j+iy 


1/2 


+ 


AT 

E 

i=0 


«i+1p-«;+1)' 


(2 


1/2 


(4.3) 


and  require  that  Eiter  be  smaller  than  an  arbitrary  tolerance  Terr.  In  all  our 
computations  we  used  10-9  <  Terr  <  10-12. 


4.1.3  Evaluation  of  u  from  U 


Once  U  has  been  evaluated  at  t  = 

(j  +  1)  At  for  x  = 

-  iAx,  i  = 

2,3,.. 

,.JV- 

3,  iV  — 

2,  a  system  of  algebraic  equations 

can  be  written  as 

-4,ri  X  TO 

b  m  —  Um 

,  m 

=  N-3 

(4.4) 

where 

c2  d2 

e2 
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0  ) 
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d-3 

e3 
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O4  64 
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d4 

e4  0 
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(4.5) 
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a-N- 4 

&AT-4 
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.4  e^r_4 

0  ... 
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0JV-3 
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CO 
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•3  djv- 3 

0  ... 
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0-N-2  b;v- 

■2  cat-2  / 

um  and  Um  are  the  vectors  containing  the  unknowns  u  and  U  at  i  =  2, ...,  N  —  2. 
Each  row  of  the  system  represents  the  finite  difference  approximation  for  the  def¬ 
inition  of  U{u)  (2.29).  a,,  bi,  ct,  d,-,  e,-  are  the  coefficients  appearing  in  front  of  u 
after  the  5-point  derivatives  are  substituted  into  (2.29),  except  for  rows  2,  3,  N-3, 
N-2,  where  these  coefficients  are  modified  to  accommodate  the  boundary  condi¬ 
tions  given  by  uxx  =  u  =  0.  After  solving  the  system  for  u  at  i  =  2, ..,  N  —  2,  we 
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use  these  conditions  (with  an  off-centered  5-point  finite-difference  approximation) 
to  obtain 


—  1f)4  —  56u^_3  +  H«_4) 


The  justification  for  using  the  boundary  conditions  uxx  =  u  =  rjx 
guarantee  that  the  mass  flux 


=  0  is  to 
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vanishes  at  the  walls,  which  can  be  verified  exactly  by  substituting  ux 
rjx  =  0  into  the  linearized  flat  bottom  momentum  equation: 

gr\x  +  u  +  ^(B -l)h2uxx 

+  ~  (^B2  —  B  — —  +  — ^  h4uxxxx  =  0 


and  obtaining  that  uxxxx  =  0,  and  therefore  that  M  =  0  at  the  boundary. 


4.1.4  Convergence  and  Stability 

No  stability  analysis  for  the  present  numerical  formulation  of  the  FN4 
model  was  done,  due  to  the  complexity  of  the  model  as  well  as  the  numerical 
scheme.  To  attain  the  desired  accuracy  in  the  model  with  relatively  fast  conver¬ 
gence,  the  Courant  number  (in  a  linear  shallow  water  theory  sense)  used  in  all 
cases  was  never  larger  than  0.3.  The  numerical  implementation  for  the  linearized 
model  proved  to  be  stable  for  all  cases  tested.  For  some  cases,  it  was  necessary  to 
filter  the  solution  as  high  frequency  oscillations  growth  appeared  near  the  points 
where  the  bottom  slope  was  discontinuous.  This  is  due  to  the  fact  that  the  FN4 
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model  contains  terms  proportional  to  high  (up  to  fourth)  derivatives  of  h  with 
respect  to  x ,  and,  for  discontinuous  bottom  slopes,  these  factors  become  singular 
and  can  introduce  spurious  high  frequency  waves  to  the  solution.  When  necessary 
we  used  a  Shapiro  (1970)  filter  with  either  8  or  16-point  average,  and  the  filter 
was  applied  every  Nf  time  steps  where  50  <  Nj  <  500,  depending  on  the  case. 

In  most  cases  where  nonlinear  effects  become  important  (very  steep,  high 
waves),  the  iteration  process  tend  to  become  slow,  or  even  diverge.  To  correct  this 
problem  we  adopted  a  relaxation  technique  in  the  iteration  process,  as  follows:  if 
with  two  iterations  the  error  tolerance  is  not  met,  we  assume  that  the  corrector 
is  overshooting  the  desired  solution  and  apply  the  formula  to  both  u  and  rj: 

St\  =  (i  -  R)f?\  +  Rfi fl  (4.7) 

where  the  relaxation  coefficient,  R  ranges  from  0  to  1.  //+1p  is  the  estimate  in 
the  previous  iteration,  and  //+1  is  the  estimate  in  the  current  iteration,  which  is 
replaced  by  the  relaxed  vector  //+1r.  The  optimal  value  of  R  strongly  depends 
the  type  of  problem.  For  most  cases,  we  used  0.2  <  R  <  1.  For  very  highly 
nonlinear  (near  breaking)  solitary  waves  it  was  necessary  to  use  R  as  small  as 
0.08  to  keep  the  solution  from  diverging.  The  number  of  iterations  necessary  for 
convergence  within  the  desired  accuracy  was  typically  less  than  6,  but  for  some 
very  near-breaking  solitary  waves  it  was  as  high  as  20. 

4.2  The  Sponge  Layer 

The  last  term  in  (4.2)  is  a  linear  friction-type  term  and  is  referred  to  as 
a  “Newtonian  cooling”  by  Israeli  and  Orszag  (1989).  Other  types  of  dissipation 
terms  such  as  viscous-type  dissipation  (analogous  to  dissipation  due.  to  viscosity 
in  the  Navier-Stokes  equations),  and  “sponge-filter”  (Israeli  and  Orszag,  1989) 
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are  also  possible,  but  it  was  found  that  the  Newtonian  cooling  was  sufficient  to 
damp  the  waves  with  efficiency  and  emulate  well  the  radiation  condition.  In 
principle,  the  only  required  rule  for  the  sponge  layer  function  /<*  is  that  it  must 
vanish  everywhere  except  near  the  boundaries,  where  the  dissipation  takes  place. 
In  practice,  a  smooth  transition  between  the  sponge  layer  and  the  interior  of  the 
domain  is  necessary  to  minimize  reflection  from  the  sponge  layer  back  into  the 
domain,  which  is  highly  undesirable.  We  choose  the  same  form  for  the  sponge 
layer  coefficient  fj(x )  as  did  Wei  (1997): 


fd(x)  = 


Sexp[(TP)N  —  ll 
exp(l)— 1 


o 


Xs  <  x  <  XL 
0  <  x  <  xs 


(4.8) 


where  xp  is  a  transformed  coordinate  defined  by 


Xp  = 


X  -  Xs 
XL  -  XS 


(4.9) 


and  5  is  a  dissipation  strength  constant.  In  all  our  computations,  N  =  2  gave 
satisfactory  results.  Reasonable  width  of  the  sponge  layer  (xp  —  xs)  and  values 
of  S  depend  on  the  wave  conditions.  Values  of  20  <  S  <  40  combined  with  the 
sponge  layer  width  of  about  three  to  four  times  the  characteristic  wavelength  gave 
satisfactory  results  for  all  cases  considered.  For  deeper  water  (shorter  waves),  the 
velocity  variable  ft  can  be  quite  small,  and  the  dissipation  term  may  not  be  able 
to  damp  the  waves  very  efficiently.  For  such  cases,  S  or  xp  —  x$  (or  both)  should 
be  increased. 


4.3  Wave  Generation  Inside  the  Domain:  Source  Function 

It  should  be  clear  that  using  sponge  layers  near  the  boundaries  rules  out 
the  possibility  of  any  type  of  wave  generation  at  the  boundary.  It  is  therefore 
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necessary  to  introduce  a  source  term  inside  the  domain  to  generate  the  desired 
waves,  which  will  be  allowed  to  propagate  towards  the  boundaries,  where  these 
waves  are  damped  by  virtue  of  the  sponge  layers.  This  is  the  role  of  the  last  term 
in  (4.1),  fs(x ,  t).  The  first  attempt  to  include  such  terms  in  Boussinesq  models  was 
made  by  Larsen  and  Dancy  (1983),  in  which  mass  is  added  and  subtracted  from 
the  domain  along  a  single  line  (point,  in  the  case  of  a  one-dimensional  model). 
Wei  (1997)  found  that  this  approach,  which  worked  well  with  the  staggered  grid 
of  Larsen  and  Dancy  (1983),  did  not  work  well  in  his  non-staggered  grid,  where 
spurious  noise  appeared  around  the  source  point.  It  was  necessary  therefore  to 
distribute  the  source  function  around  a  certain  neighborhood  of  the  source.  In 
the  present  formulation,  we  closely  follow  the  approach  of  Wei  (1997),  in  which 
the  source  function  is  assumed  to  be  distributed  as  a  Gaussian  shape,  making  the 
appropriate  modifications  to  account  for  the  added  complexity  of  the  model.  The 
formulation  for  the  source  function  presented  next  is  one-dimensional,  but  can  be 
extended  to  two  dimensions  is  a  straightforward  manner. 

If  the  local  water  depth  at  the  source  region  is  constant,  h,  and  we  want 
to  generate  regular  waves  with  angular  frequency  uj,  the  source  function  can  be 
written  as: 

fa(x ,  t)  =  Ds  exp[— /?*(#  -  xs)2]  sin (wf),  (4-10) 


where  xs  is  the  center  of  source  function,  [3S  determines  how  focused  the  source 
function  is,  and  Ds  is  the  magnitude  of  the  source  function.  Assuming  that  the 
generated  wave  have  small  amplitude,  we  can  use  the  linearized  version  of  the  FN4 
model  and  derive  an  exact  analytical  expression  for  Ds  by  using  Green’s  function 
theory  (see  Appendix  C),  to  obtain: 


Ds 


_ Vo _ 

uasI\  [l  +  C3(kh)2  +  C4(kh)4Y 


(4.11) 
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where  k  is  the  wave  number  computed  from  the  linear  dispersion  relationship  for 
the  present  model  (see  chapter  3),  of  the  wave  generated  at  the  source  function 
when  it  is  away  from  the  source  region,  C3  and  C4  are  constants  given  in  chapter 
3,  as  and  I\  are  given  in  appendix  C. 

Although  the  gaussian  shape  parameter  j3s  is  arbitrary,  in  practice  its  value 
has  great  influence  on  how  well  the  source  function  can  generate  the  desired  waves. 
Ideally,  / 3  should  be  as  large  as  possible,  so  that  the  source  function  would  be  more 
localized.  However,  it  turns  out  that  if  the  source  region  is  too  narrow  (large  /3S), 
the  waves  generated  can  be  quite  distorted  and  noise  may  also  appear  when  the 
waves  not  are  small  amplitude  (see  next  subsection).  Defining  the  width  of  the 
source  region  Ws  to  be  the  distance  between  two  coordinates  (equidistant  from 
the  source  center)  where  exp[— {3s(x  —  xs)2]  is  equal  to  e~5,  we  can  write: 

W.  =  2y/b/P..  (4.12) 

By  trial  and  error,  it  was  found  that,  for  regular  waves,  a  source  with  width 
Ws  approximately  equal  to  the  wave  length,  gives  satisfactory  results  for  waves 
within  a  wide  range  of  amplitudes  and  wavenumber.  Sensitivity  tests  for  Ws  were 
performed  and  are  shown  next. 

4.3.1  Tests 

Tests  of  the  source  function  and  sponge  layers  are  now  presented.  Figure 

4.1  shows  schematically  the  location  of  the  source  and  the  sponge  layers.  In  all 
the  cases  a  wave  with  period  T  =  1.5s,  and  wavelength  L  =  3.35m  was  used.  The 
water  depth  is  h  =  1  m,  sponge  layers  widths  are  3 L,  and  sponge  layer  strength 
S  =  30.  Six  cases  were  considered:  amplitudes  do  =  0.01m,  a 0  =  0.05m,  and 
a0  =  0.1m,  each  for  source  region  with  Ws  =  0.25L  and  Ws  =  L.  Figure  4.2 
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Figure  4.1:  Sketch  of  the  domain  for  source  function  and  sponge  layers  tests. 
Domain  length  Lx  =  80m  was  used  throughout. 

shows  results  at  various  times  prior  to  “steady-state”,  for  ao  =  0.01m,  and  both 
source  widths.  Notice  that  for  this  particular  case  of  very  small  amplitude  waves, 
the  source  function  worked  quite  well  in  generating  the  desired  waves,  regardless 
of  the  width  Ws.  With  ao  =  0.05m,  a  more  “nonlinear  wave”,  4.3  shows  that 
the  narrower  source  region  was  not  able  to  generate  clean  waves,  and  it  can  be 
seen  that  higher  frequencies  contaminate  the  solution.  Notice  that  these  higher 
frequencies  are  not  necessarily  near  the  Nyquist  (instability-type)  frequency.  It 
seems  that  the  higher  harmonics  generated  by  the  model  are  not  behaving  as  we 
desire  as  they  are  generated  in  the  source  function  (which  is  derived  from  linear 
theory).  Using  Ws  =  L,  however,  this  problem  is  corrected  for  this  amplitude. 
It  seems  that  a  good  resolution  of  the  source  function,  that  is  a  larger  Ws,  is 
necessary  to  keep  undesirable  higher  frequencies  from  appearing  in  the  solution. 
4.3  shows  similar  results  for  even  higher  waves  a0  =  0.1m.  With  Ws  =  0.25L,  the 
solution  became  unstable  right  after  t  =  4T,  whereas  for  Ws  =  L,  the  resulting 
waves  are  still  quite  clean. 
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Figure  4.2:  Waves  with  a0  =  0.01m  generated  by  a  source  function  at  x  =  40m 
with  Ws  =  0.25L  (upper  panel),  and  Ws  =  L  (lower  panel),  h  =  lm, 
T  =  1.5s,  xl  —  xs  =  3 X,  S  =  30. 


Figure  4.3:  Waves  with  a0  =  0.05m  generated  by  a  source  function  at  x  =  40m 
with  Ws  =  0.25 L  (upper  panel),  and  W3  =  L  (lower  panel),  h  =  1  m, 
T  =  1.5s,  xl  —  xs  =  3 L,  S  =  30. 


Figure  4.4:  Waves  with  ao  =  0.1m  generated  by  a  source  function  at  x  =  40m 
with  Ws  =  0.25 L  (upper  panel),  and  Ws  =  L  (lower  panel),  h  =  lm, 
T  =  1.5s,  xl  —  xs  =  3L,  S  =  30. 


Chapter  5 


THE  SOLITARY  WAVE 


The  phenomenon  known  as  the  solitary  wave  consists  of  a  limiting  wave 
form  with  a  single  crest  which  propagates  in  fairly  shallow  water  of  constant  depth, 
and  where  the  nonlinear  and  dispersive  effects  counterbalance  each  other  yielding 
a  permanent  form  solution.  In  this  chapter  we  study  solitary  waves  solutions  of 
the  FN4  model,  and  compare  it  to  other  models  including  extremely  accurate 
solutions  of  the  full  boundary  value  problem  (Tanaka,  1986),  Green  and  Naghdi 
(1976)  (GN)  type  models,  and  the  WKGS  model. 

Many  authors  have  found  approximate  solutions  for  the  solitary  wave,  in¬ 
cluding  the  early  works  of  Boussinesq  (1871)  and  Korteweg  and  deVries  (1895). 
Fenton  (1972)  developed  a  model  based  on  a  perturbation  expansion  around  the 
basic  shallow  water  wave  theory.  His  expansion  includes  terms  up  9th  order  and, 
at  the  first  three  orders,  recover  the  models  of  Boussinesq  (1871),  Laitone  (1960), 
and  Grimshaw  (1971).  Longuet-Higgins  and  Fenton  (1974)  used  conservation  of 
integral  quantities  such  as  mass  and  energy  to  arrive  at  extremely  accurate  re¬ 
lationships  between  several  solitary  wave  properties,  such  as  the  wave  height, 
energy,  mass,  wave  Froude  number  Fr  (nondimensional  wave  speed),  etc.  They 
also  proved  that  the  solitaxy  wave  with  maximum  wave  height  does  not  corre¬ 
spond  to  the  one  with  maximum  fluid  velocity  at  the  crest,  or  maximum  mass. 
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More  recently,  in  a  study  of  the  stability  of  solitary  waves,  Tanaka  (1986)  devel¬ 
oped  an  accurate  solution  scheme  to  the  full  boundary  value  problem  for  solitary 
waves.  Throughout  this  chapter  we  will  use  this  solution  as  the  “exact”  solution 
in  our  comparisons.  Shields  and  Webster  (1988)  studied  the  accuracy  of  solitary 
wave  properties  of  the  first  3  levels  of  the  GN  models  (referred  hereafter  as  GN1, 
GN2,  and  GN3).  An  nth  level  GN  model  approximates  the  horizontal  velocity 
by  an  ( n  —  l)th  polynomial,  and  the  vertical  velocity  by  an  nth  order  polynomial. 
GN1  recovers  the  model  by  Serre  (1953),  as  shown  by  Kirby  (1997).  Shields  and 
Webster  (1988)  derived  a  GN2  set  of  equations  for  unsteady  flow  over  an  uneven 
bottom,  and  a  GN3  model,  for  one-dimensional  steady  flow  over  a  flat  bottom. 

In  the  next  sections,  we  compare  several  properties  of  the  solitary  wave 
solution  of  the  FN4  model  with  the  exact  solution,  numerical  WKGS  solutions, 
and  also  with  GN  solutions  given  by  Shields  and  Webster  (1988). 

5.1  Linear  Analytical  Asymptotic  Solution 

At  the  tail  of  the  solitary  wave  (away  from  the  crest)  the  free  surface 
elevation  r/  is  very  small,  and  we  expect  that  the  linearized  set  of  equations  should 
describe  the  shape  of  the  wave  with  good  accuracy.  In  a  reference  frame  moving 
with  the  wave  at  nondimensional  wave  at  speed  Fr  =  c/y/gK ,  we  can  write  the 
following  boundary  value  problem  for  the  wave  field  far  from  the  crest  (located  at 
x  =  0)  in  (x,  z): 


v24>  = 

0 

(5.1) 

(px  - 

— Fr  x  — y  00 

(5.2) 

< pz  = 

0  z-- 1 

(5.3) 

<pz  = 

—F2(pxx  Z  =  0 

(5.4) 
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The  solution  to  the  system  above  is 


4>  =  K\e <i~ix  cos  27  z  —  Frx 
Substituting  (5.5)  into  (5.4)  gives 

tan  27  2 

~^r=F- 


(5.5) 


(5.6) 


The  exact  solution  for  the  free  surface  elevation  7  far  from  the  crest  is  of  the  form 

77  =  JT2e2^  (5.7) 

The  parameter  7  is  referred  to  as  the  straining  parameter,  and  (5.6)  is  directly 
related  to  the  exact  dispersion  relationship  in  linear  wave  theory. 

For  the  present  model,  the  2  equations  corresponding  to  the  system  (5.1- 
5.4),  in  terms  of  the  modified  velocity  variable  u  are: 


7  —  F r  fX  CqUxx  C4  'U’xxxx^J 

F\ rTj  =  (u  /2  ClUxx  +  ^'2^'xxxo^j 


(5.8) 

(5.9) 


where  Ci,  C2,  C3,  C4  are  defined  in  Appendix  B.  We  now  assume  the  solution 

(5.10) 


u  =  e 


—Ax 


Substituting  (5.10)  into  (5. 8, 5. 9),  solving  for  /FX2  and  keeping  the  relevant  root, 
we  obtain 


A*2  A2  = 


(C3Fr2  -  Cr) 

1  - 

(C,F.2  -  C.) 

1 2  +  4  | 

(a 

-  F,!)  (c,f!  -  c2) 

1/2 

2  (C4Fr2  -  ( 

\ 

(5.11) 
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To  recover  Nwogu  (1993),  refer  to  chapter  3.  The  straining  parameter  7  is  related 
to  A  as  27  =  —A. 

The  expressions  for  the  relation  between  the  straining  parameter  and  the 
Froude  number  for  GN2,  and  GN3  (Shields,  1986  and  Shields  and  Webster,  1988) 
are: 

-  12F"2  -  4^/9 F-4  -  33 F-2  +  124^  (5.12) 

and 

16Fr276  +  60  (l  -  9F2)  74  -  20  (39  -  144F2)  72  +  1575  (l  -  F2)  =  0  (5.13) 

respectively. 

Figure  5.1  shows  comparison  between  the  percentage  error  to  the  exact 
solution  of  the  transcendental  equation  (5.6)  for  7,  between  the  present  model, 
Nwogu’s  model,  GN2,  and  GN3.  Although  all  models  have  relatively  small  errors, 
the  present  model  is  much  more  accurate  than  all  the  others  by  at  least  an  order 
of  magnitude.  Notice  that  compared  to  GN2  and  Nwogu’s  model,  the  difference 
in  the  errors  is  of  at  least  five  order  of  magnitude. 

5.2  Solitary  Waves  with  Permanent  Form 

In  this  section,  the  numerical  scheme  described  in  chapter  4  is  used  to  com¬ 
pute  several  approximate  solitary  wave  solutions  to  the  fully  nonlinear  models  FN4 
and  WKGS.  The  initial  condition  used  for  the  model  was  constructed  from  the 
computer  program  by  Tanaka  (1986)  in  the  following  manner:  for  the  smallest 
computed  wave  with  amplitude  rjmax  «  0.2,  r]  and  u  ( ua  in  the  case  of  WKGS) 
were  obtained  from  Tanaka’s  exact  solution  and  used  as  initial  condition  for  FN4 
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Figure  5.1:  Percentage  error  to  exact  solution  in  straining  parameter.  Present 
model  (full),  Nwogu  (dash),  GN2  (dot),  GN3  (dash-dot) 
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and  WKGS.  After  the  solution  reached  permanent  form,  it  was  multiplied  by  a 
factor  slightly  larger  than  one  (typically  1.05)  and  this  re-scaled  wave  was  used 
as  the  initial  condition  for  the  next  case.  This  procedure  was  repeated  until  the 
desired  range  of  solitary  waves  was  covered,  and  proved  to  be  more  efficient  than 
using  Tanaka’s  solution  as  the  initial  condition  for  all  amplitudes.  Since  Tanaka’s 
and  each  consecutive  re-scaled  initial  condition  do  not  satisfy  the  approximate 
equations,  there  is  a  transient  period  while  the  solution  is  not  a  permanent  form 
solitary  wave,  but  has  a  dispersive  tail  of  shorter  waves  left  behind.  Since  these 
shorter  waves  travel  with  phase  speeds  which  are  smaller  than  the  “main  wave”, 
eventually  the  tail  is  left  far  behind  and  does  not  interfere  with  the  solitary  wave, 
which,  at  this  point,  can  propagate  with  permanent  form.  The  time  required  for 
the  solution  to  achieve  a  permanent  form  solitary  wave  depends  on  the  initial  con¬ 
dition.  High  amplitude  initial  conditions  will  reach  permanent  form  more  quickly, 
since  the  primary  wave  moves  much  faster  than  the  tail.  Smaller  amplitude  waves 
will  have  less  amplitude  dispersion  and  it  will  take  longer  for  the  permanent  form 
solitary  wave  to  separate  from  the  tail.  After  each  solution  reaches  permanent 
form  it  is  straightforward  to  obtain  properties  such  as  the  Froude  number,  velocity 
profiles,  mass,  energy,  etc.  No  filtering  was  necessary  during  these  computations, 
although  for  higher  waves  the  under-relaxation  parameter  had  to  be  as  small  as 
r  =  0.08  for  the  solution  to  converge  with  error  tolerance  in  the  iteration  typi¬ 
cally  10~12  <  Terr  <  10-9.  For  the  permanent  form  solitary  waves  we  used  grid 
spacing  Ax  =  0.1  h  for  waves  with  amplitude  0  <  rjmax  <  0.4/i,  Ax  =  0.05h  for 
OAh  <  rjmax  <  0.7 h,  and  Ax  =  0.05/i  for  rimax  >  0.7 h,  where  h  =  1.  We  used  At 
such  that  the  Courant  number  was  always  below  0.2  (for  accuracy  purposes). 


We  now  present  the  following  nondimensional  quantities  for  the  FN4  model, 
where  the  scales  for  the  basic  variables  are  x  =  x'h,  2  =  z'h ,  t  =  t' f \JgJh , 
u  =  u'\fgh.  and  the  primes  denote  nondimensional  quantities.  The  primes  in  the 
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formulas  below  are  dropped  for  the  sake  of  notation  clarity.  The  total  mass  of  the 
solitary  wave  above  the  still  water  level  is  given  by: 


M  = 


(5.14) 


The  potential  energy  is: 


/+°°  1  „ 
-V2dx 

-oo  Z 


The  kinetic  energy  is: 


K  =  \C  j” 


(5.15) 


(5.16) 


where  £  =  (1  +  z),  H  =  1  +  77, 

u(C)  =  u  +  i/i2  (5  -  C2)  uxx  +  ^fi4  -  C2  -  y  +  ux 


(5.17) 


and 


M>«)  =  +  /'‘1  +  )(’)  5. 


(5.18) 


After  substitution  of  (5.17)  and  (5.18)  into  (5.16)  and  retaining  terms  of  up  to 
0(/i4),  we  obtain 

K  =  ^  J  {/i2  [u2  +  [B  -  C2)  uxuxx\ 

+  [(B  -  C2)  B  -  ^  (l>  -  C4)]  Uxuxxxx 

+  4  Mxir] 

+  (C^)2  +  c2  (b  -  ^c2)  uxuxxx}  d(dx.  (5.19) 


The  integral  in  £  can  be  evaluated  analytically,  and  the  integrals  in  x  are  computed 
using  Bode’s  rule,  which  is  accurate  to  0( Ax6),  within  a  range  containing  the 
solitary  wave  and  where  the  free  surface  at  the  extremes  is  negligibly  small. 
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We  define  the  quantity 

us  =  l  -  (tic  -  Fr)2,  (5.20) 

where  uc  is  the  particle  velocity  at  the  crest,  computed  from  (5.17)  by  locating  the 
position  of  the  crest  in  x  and  computing  u((  =  H )  at  that  x  location.  As  the  wave 
amplitude  varies  from  0  to  its  limiting  value,  in  which  uc  =  Fr,  the  parameter  us 
goes  from  0  to  1. 

The  speed  of  each  wave  was  computed  by  letting  an  already  permanent  form 
solution  propagate  over  a  distance  of  around  500  times  the  water  depth,  recording 
the  difference  between  the  crest  location  xc  before  and  after  this  interval  dt  and 
computing 


The  exact  location  of  the  wave  crest  could  not  be  obtained  directly  from 
the  computations,  since  only  by  virtue  of  luck  the  crest  was  located  exactly  at  one 
of  the  grid  points.  The  location  of  the  crest  was  determined  by  fitting  a  4th  order 
polynomial  to  the  free  surface  around  the  crest.  The  peak  value  and  x  location 
were  then  computed  using  the  fitted  polynomial.  We  used  this  same  approach  to 
compute  u  and  its  x  derivatives  at  and  underneath  the  crest. 

Figure  5.2  shows  computations  of  the  free  surface  elevation  of  half  of  a  solitary 
with  Fr  =  1.266  for  FN4,  WKGS,  GN1,  GN2,  GN3,  and  the  exact  solution.  The 
three  GN  models  are  plotted  with  dotted  lines,  with  GN1  and  GN2  . marked  with 
labels.  Notice  that,  of  all  models,  GN3  has  the  best  match  with  the  exact  solution. 
FN4  is  also  fairly  close  to  the  exact  solution,  but  WKGS  strongly  overpredicts  the 
wave  height,  and  slightly  underpredicts  the  tail.  GN  models  tend  to  underpredict 
the  height  and  overpredict  the  tail.  Figure  5.3  shows  the  same  model  comparisons 
as  in  Figure  5.2,  except  for  GN1  and  GN2  whose  solutions  were  not  available.  In 
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this  case  the  maximum  wave  height  is  kept  constant  for  all  models.  Again,  GN3 
has  the  best  shape  compared  to  the  exact  solution.  WKGS  solution  compares 
better  with  the  exact  solution  than  in  the  previous  case  ( Fr  kept  constant),  but 
FN4  still  compares  better  with  the  exact  solution  than  does  WKGS.  Notice  also 
the  difference  in  wave  speed  Fr  for  each  model.  In  Figure  5.4  the  vertical  profiles 
of  the  horizontal  velocity  are  shown  for  the  exact  solution,  FN4,  and  WKGS,  for 
the  waves  shown  in  Figure  5.2,  and  it  can  be  seen  that  the  0(/i4)  model  has  a 
more  accurate  kinematics  representation  than  the  0(fi2)  model,  confirming  what 
we  have  already  shown  for  linear  theory.  Unfortunately  it  was  not  possible  to 
obtain  GN  vertical  profiles,  but  it  can  be  speculated  that  it  would  not  be  able 
to  predict  this  property  as  accurately  as  the  FN4,  since  it  assumes  the  horizontal 
velocity  to  be  only  a  2nd  order  polynomial. 

Figure  5.5  shows  the  relationship  between  the  speed  and  amplitude  of  a 
wide  range  of  solitary  waves  for  several  models.  Notice  that  once  again  GN3  has 
the  closest  solution  to  the  exact  one.  FN4  slightly  underpredicts  the  wave  speed 
for  a  given  amplitude,  whereas  the  deviation  in  WKGS  is  of  an  order  of  magnitude 
higher.  GN2  and  especially  GN1  overpredict  the  wave  speed  throughout  the  range 
tested.  It  is  important  to  keep  in  mind  that  as  the  wave  approach  the  limiting 
value,  the  crest  becomes  extremely  sharp  (with  the  limiting  wave  having  a  crest 
forming  an  angle  of  120°),  which  makes  it  difficult  for  the  finite  difference  scheme 
of  FN4  to  resolve  the  wave  well  near  the  crest,  since  the  model  has  up  to  5th  order 
derivatives  in  x. 

Figure  5.6  shows  computations  of  the  parameter  uis  as  a  function  of  the 
wave  speed.  In  this  case,  the  FN4  model  is  the  closest  to  the  exact  solution.  This 
is  not  surprising  if  one  recalls  that  us  is  directly  related  to  the  horizontal  fluid 
velocity  at  the  crest,  and  that  the  FN4  has  a  4th  order  polynomial  representation 
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Figure  5.2:  Shape  of  solitary  waves  with  fixed  Fr  =  1.266.  Exact  (full),  FN4 
(dash),  WKGS  (dash-dot),  GN1,  GN2,  GN3  (dot) 
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Figure  5.3:  Shape  of  solitary  waves  with  fixed  amplitude  r/max  =  0.65.  Exact: 

Fr  =  1.265  (full),  FN4:  Fr  =  1.262  (dash),  WKGS:  Fr  =  1.245 
(dash-dot),  GN3:  Fr  =  1.266  (dot) 


70 


Figure  5.5:  Solitary  wave  amplitude  vs.  phase  speed.  Exact  (full),  FN4  (dash), 
WKGS  (dash-dot),  GN1,  GN2,  GN3  (dot) 

of  the  vertical  profile  of  the  horizontal  velocity,  whereas,  as  already  observed,  only 
a  2nd  order  polynomial  is  assumed  in  both  GN3  and  WKGS.  In  the  next  figures, 
GN  solutions  were  not  available.  Figures  5.7,  5.8,  and  5.9  show  plots  of  the  mass, 
kinetic  energy,  and  potential  energy  of  solitary  waves  against  the  wave  speed,  for 
the  exact  solution,  FN4,  and  WKGS.  All  three  properties  show  a  similar  behavior 
to  the  wave  amplitude  (Figure  5.5)  when  plotted  against  Fr.  In  Figure  5.7,  WKGS 
agrees  with  the  exact  solution  better  than  FN4,  but  this  is  only  a  coincidence, 
as  the  overprediction  of  the  wave  crest  counterbalance  the  underprediction  of  the 
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Figure  5.7:  Mass  vs.  phase  speed  for  solitary  waves.  Exact  (full),  FN4  (dash), 
WKGS  (dash-dot) 

wave  tail.  A  similar  effect  happens  with  the  kinetic  energy  (Figure  5.8),  where 
the  plots  of  the  2  models  coincidently  are  on  top  of  each  other.  The  potential 
energy  (Figure  5.9)  calculations  for  model  FN4  has  better  agreement  to  the  exact 
solution  than  it  has  for  the  WKGS  model,  which  confirms  that  the  good  agreement 
of  WKGS  in  Figure  5.5  was  by  virtue  of  luck,  since  both  the  mass  and  the  potential 
energy  are  only  dependent  of  the  free  surface  elevation. 
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5.2.1  Discussion 


From  Figures  5.1,  5.2,  5.3,  5.5,  5.6,  it  is  clear  that  the  FN4  model  has  a 
better  asymptote  (linear)  agreement  with  the  exact  solution  than  does  GN3,  but, 
with  the  exception  of  the  parameter  us  (related  to  the  velocity  at  the  crest  of  the 
wave),  in  all  other  nonlinear  properties,  GN3  has  a  better  agreement  than  FN4. 
This  may  seem  somewhat  surprising  since  GN3  approximates  the  horizontal  veloc¬ 
ity  by  a  second  order  polynomial  (two  orders  lower  than  the  FN4  model)  and  the 
vertical  velocity  by  a  third  order  polynomial  (same  as  the  FN4  model),  and  a  more 
careful  study  is  needed  to  explain  these  discrepancies.  Nevertheless,  we  make  the 
following  conjectures:  FN4  satisfies  mass  conservation  and  all  boundary  condi¬ 
tions  in  an  approximate  sense,  consistent  with  the  level  of  approximation  of  the 
velocity  field.  GN3  satisfies  mass  conservation  and  the  kinematic  boundary  con¬ 
ditions  exactly.  The  coefficients  for  the  velocity  variable  in  FN4  are  derived  such 
that  the  linear  dispersion  relationship  is  extremely  accurate,  and  no  optimization 
is  done  considering  that  the  free  surface  displacement  is  finite  or  including  nonlin¬ 
ear  terms.  The  advantage  of  this  approach  is  that  the  model  is  simple  in  the  sense 
that  there  is  only  one  dependent  variable  describing  the  internal  kinematics.  For 
small  amplitude  waves  the  FN4  model  is  capable  of  extremely  accurate  results  for 
a  wide  range  of  water  depths,  but  as  the  nonlinear  terms  become  more  important, 
and  the  actual  free  surface  deviates  considerably  from  the  still  water  level  (for 
which  the  model  was  derived  to  perform  its  best),  errors  start  to  increase.  On  the 
other  hand,  GN3  is  derived  without  any  assumption,  and  no  optimization  is  made 
“a  priori”.  Rather,  the  coefficients  for  the  velocity  polynomial  are  derived  by  a 
minimization  of  the  errors  in  the  momentum  equation  over  the  entire  actual  water 
depth,  which  is  changing  with  time.  As  a  consequence,  GN3  has  three  dependent 
variables  describing  each  component  of  the  velocity  field,  besides  the  free  surface 
elevation.  For  a  2-dimensional  problem,  GN3  would  have  seven  coupled  evolution 
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Figure  5.10:  Sketch  of  shoaling  solitary  wave, 
equations,  whereas  FN4  would  have  three. 

5.3  Shoaling  Solitary  Wave 

In  this  section  we  use  permanent  form  solitary  wave  solutions  for  the  FN4 
and  WKGS  models  as  initial  conditions  for  a  more  complex  domain  which  includes 
a  slope  where  the  waves  shoal  and  eventually  break.  The  two  Boussinesq  models 
are  compared  to  the  Boundary  Element  Method  (BEM)  results  using  Grilli  et  al. 
(1989)  formulation,  which  is  a  solution  to  the  full  boundary  value  problem,  and 
can  be  regarded  as  an  “exact  solution”  for  the  sake  of  comparison  with  Boussinesq 
models.  Comparisons  between  WKGS,  Nwogu  (1993),  and  BEM  for  several  cases 
can  be  found  in  Wei  et  al.  (1995).  Here,  we  concentrate  in  a  single  case  of  a  wave 
with  nondimensional  amplitude  r)max/ho  =  0.2,  where  h0  is  the  water  depth  before 
the  slope  s  =  1/35.  A  sketch  of  the  problem  is  illustrated  in  Figure  5.10. 

Figure  5.11a  shows  the  free  surface  elevation  at  4  different  times  as  the  wave 
shoals  on  the  slope.  The  nondimensional  x'  =  x/ho  has  its  origin  at  the  toe  of  the 
slope,  and  the  nondimensional  time  t1  =  t(halg)~ll 2  has  its  origin  at  the  instant 
the  wave  crest  passes  through  x1  =  0.  This  definition  was  used  to  synchronize  the 
3  models’  results.  Although  the  differences  are  subtle,  it  is  possible  to  see  that 


78 


(a) 


Figure  5.11:  (a)  Solitary  wave  shape  at  t'  =  t\  =  39.98,  ^2  =  53.19,  =  61.13, 

t4  =  66.89.  (b)  Crest  speed  (top  curves)  and  fluid  velocity  (bottom 
curves).  Circle  denotes  breaking  point.  BEM  (full),  FN4  (dash), 
WKGS  (dash-dot) 

the  FN4  wave  is  closer  to  the  BEM  than  the  WKGS  solution,  in  particular  near 
the  wave  crest  at  t3  and  t4  (shown  in  finer  detail  in  Figure  5.12).  Figure  5.11b 
shows  computations  of  the  crest  speed  and  the  fluid  velocity  at  the  crest.  The 
location  of  breaking,  defined  as  the  point  where  a  vertical  tangent  is  formed  at  the 
face  of  the  wave  from  the  BEM,  is  also  shown.  The  crest  speed  was  computed  by 
locating  and  recording  the  wave  crest  position  xc  at  every  time  step  t,  in  a  similar 
manner  to  what  was  done  in  the  previous  section,  then  computing  the  speed  by 

Co  =  ^  (5.22) 
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Figure  5.12:  (a)  Solitary  wave  shape  at  t'  =  t4  =  66.89.  BEM  (full),  FN4  (dash), 
WKGS  (dash-dot) 
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Near  breaking,  the  computations  of  Cc  and  Vc  were  smoothed  with  a  12th  order 
polynomial,  since  the  estimations  given  by  (5.22)  exhibited  some  jitter.  The  model 
ran  only  up  to  the  point  where  Cc  =  Vc ,  and  the  values  beyond  that  point  on 
Figure  5.11b  are  due  to  the  polynomial  extrapolation.  The  fluid  velocity  at  the 
crest  which,  due  to  the  presence  of  the  slope,  has  both  horizontal  and  vertical 
components,  is  given  by: 

Ve  =  \[u?c  +  (5-23) 

where  uc  and  wc  can  be  evaluated  from  the  formula  for  the  velocity  profiles  given 
below: 

w  =  u  +  H2  [{Ah  —  C)/21ar  +  2 hx(Ah  —  C)/22  +  ( Bh 2  —  C2)/22a;] 

+  /i4  [{Ah  -  0/41.  +  2 hx(Ah  -  0/42  +  ( Bh 2  -  CO/42, 

+  {Ah  -  0/43,  +  2 hx{Ah  -  0/44  +  ( Bh 2  -  e)Uix 
+  3hx{Bh2-(2)f45  +  {Ch3-C3)f45x 

+  4 hx{Ch3  —  C3)/i6  +  {Dh4  —  C4)/46;r]  j  (5.24) 

w  =  —  [n2  (/21  +  2^/22)  +  (/41  +  2C/42 

+  /43  +  2^/44  +  3O/45  +  4(3/46)]  •  (5.25) 

where,  again,  nondimensional  quantities  are  implied  but  primes  have  been  omit¬ 
ted  for  notation  simplicity.  Model  FN4  performs  better  than  WKGS  in  predicting 
both  the  phase  speed  and  the  crest  fluid  velocity.  Finally,  vertical  profiles  of  the 
horizontal  velocity  under  the  crest  are  shown  in  Figure  5.13  at  three  different  lo¬ 
cations,  including  one  near  the  breaking  point  as  predicted  by  BEM.  As  expected, 
the  FN4  model  performs  better  than  the  WKGS  due  to  its  better  representation 
of  the  internal  kinematics. 
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Chapter  6 


COMPARISONS  WITH  LABORATORY 
MEASUREMENTS 


It  is  well  known  that  regular  waves  decompose  into  higher  frequency  free 
waves  as  they  propagate  past  a  submerged  bar,  as  shown  in  experimental  work 
by  Beji  and  Battjes  (1993),  Luth  et  al.  (1994),  and  Ohyama  et  al.  (1994).  The 
basic  mechanism  is  as  follows:  as  the  waves  propagate  onto  the  front  slope  of 
the  bar,  nonlinear  interactions  transfer  energy  from  the  leading  wave  component 
(primary  wave)  to  higher  harmonics,  causing  the  wave  to  become  steeper  and 
also  asymmetric  (pitched  forward).  After  the  peak  of  the  bar  is  reached  (say 
no  breaking  occurs),  and  the  bottom  slope  becomes  negative  (depth  increases), 
the  nonlinear  coupling  (forcing)  of  the  higher  harmonics  with  the  fundamental 
wave  becomes  progressively  weaker,  and,  from  higher  to  lower  harmonics,  each  of 
the  Fourier  components  are  released  as  free  waves  with  their  own  bound  higher 
harmonics.  Of  course,  since  the  waves  after  the  bar  travel  with  different  speeds, 
the  process  can  be  fairly  complicated  with  some  waves  overtaking  others,  and 
involving  nonlinear  interactions.  It  is  clear  therefore  that  wave  propagation  over 
a  submerged  bar  is  a  quite  demanding  test  for  Boussinesq-type  models,  as  it 
requires  that  the  model  predict  the  nonlinear  harmonic  generation  well,  and  also 
that  the  released  shorter  waves  (behind  the  bar)  have  an  accurate  speed,  which 
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may  not  happen  even  if  the  model  predicts  well  the  speed  of  the  primary  waves 
before  they  reach  the  bar. 

Comparisons  between  several  weakly  nonlinear  Boussinesq-type  models  and 
experimental  data  by  Beji  and  Battjes  (1993)  and  Luth  et  al.  (1994)  of  waves  prop¬ 
agating  over  a  submerged  bar  were  presented  by  Dingemans  (1994).  In  general, 
the  models  performed  relatively  well  for  the  longer,  lower  amplitude  waves,  but 
all  were  fairly  inaccurate  for  the  shorter,  more  nonlinear  waves,  especially  behind 
the  bar. 

Comparisons  between  the  extended  Boussinesq  model  by  Nwogu  (1993), 
among  other  types  of  models,  and  experimental  data,  also  for  waves  passing  over 
a  submerged  bar,  were  presented  by  Ohyama  et  al.  (1994),  and  the  results  were 
similar  to  the  comparisons  made  by  Dingemans  (1994),  that  is,  the  model  poorly 
predicted  waves  behind  the  shoal  for  the  shorter,  higher  wave  cases. 

In  this  chapter  we  compare  the  FN4  model  with  three  laboratory  exper¬ 
imental  data  sets  of  regular  waves  propagating  in  a  one-dimensional  wave  flume 
and  over  a  submerged  bar:  Beji  and  Battjes  (1993),  Luth  et  al.  (1994),  and 
Ohyama  et  al.  (1994).  We  also  show  comparisons  of  the  WKGS  model  with  the 
same  data  sets.  The  models’  comparisons  with  the  data  are  done  in  three  different 
manners:  plots  of  free  surface  time  series  at  fixed  locations,  spatial  plots  of  Fourier 
components  of  the  time  series,  and  a  quantitative  estimation  of  accuracy  defined 
by: 

£  [y(j)  -  yd{j )]2 

di  =  1  -  n2  J=ni -  (6.1) 

£  [| y(j)  -  yd\  + 1 yd{j)  -  yd\f 

j~n  i 

where  d,-  is  an  index  of  agreement  proposed  by  Wilmott  (1981)  for  the  ith  wave 
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gauge,  and  where  ni  and  n2  cover  a  full  wave  period  in  the  time  series. 

Ud{j)  are  the  measured  data  to  be  compared  with,  y(j)  are  the  predicted 
values  from  the  model,  and  yd  is  the  mean  value  of  yd{j)-  A  perfect  agreement 
between  data  and  model  corresponds  to  d t-  =  1,  while  a  complete  disagreement 
results  in  d,  =  0.  In  all  the  numerical  simulations  we  used  Ax  =  0.025m  and  kept 
the  Courant  number  below  0.3.  The  sponge  layer  strength  and  width  used  were 
S  =  30  and  xl  —  xs  =  3 L,  respectively.  The  width  of  the  source  function  used 
was  Ws  =  L,  where  L  is  the  incident  wave  length.  The  details  of  the  experiments 
and  comparisons  with  the  models  are  presented  in  the  sections  below. 

6.1  The  Delft  Hydraulics  Experiments 

The  experiments  performed  by  Beji  and  Battjes  (1993)  and  Luth  et  al. 
(1994)  have  the  same  geometric  characteristics,  except  for  the  length  scale  in 
Luth  et  al.  (1994),  which  is  twice  as  large  as  in  Beji  and  Battjes  (1993).  In  Luth 
et  al.  (1994)  all  gauge  locations  used  in  Beji  and  Battjes  (1993)  were  repeated,  and 
another  run  of  measurements  was  performed  with  the  gauges  at  different  locations. 
For  the  sake  of  consistency  with  the  study  by  Dingemans  (1994),  we  re-scale  all 
measurements  to  the  scales  used  in  Beji  and  Battjes  (1993).  The  layout  of  the 
experimental  set-up  with  the  locations  of  the  measurement  stations  (to  which  we 
refer  by  their  location,  e.g.  gauge  2.0m,  gauge  15.7m,  etc)  and  the  geometry  of 
the  flume  are  illustrated  in  Figure  6.1.  In  the  present  work  we  use  the  data  from 
Luth  et  al.  (1994),  since  in  that  experiment  active  wave  absorption  was  used  at 
the  end  of  the  flume  and  both  reflection  and  bound  long  waves  were  monitored 
during  the  experiment. 

Three  sets  of  data  were  collected  using  different  incident  wave  conditions. 
We  refer  to  these  data  sets  as  cases  (a),  (b),  and  (c).  In  case  (b),  wave  breaking 
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Gage  locations 

2.0  4.0  5.7  10.5  12.5  13.5  14.5  15.7  17.3  19.0  21.0  23.0 


Figure  6.1:  Sketch  of  wave  flume  of  Delft  experiments.  All  dimensions  in  (m) 


Table  6.1:  Incident  wave  characteristics  for  the  Delft  experiments. 


Case  (a) 

Case  (c) 

Wave  amplitude  (m) 

0.01 

0.0205 

Wave  period  ( s ) 

2.02 

1.01 

O 

O 

.-Sd 

III 

=1 

0.67 

1.69 

o 

o 

<3 

III 

0.025 

0.051 

occurred  on  the  crest  of  the  shoal,  and  therefore  these  data  were  disregarded,  since 
the  present  model  does  not  include  any  breaking  mechanism.  The  incident  wave 
characteristics  for  cases  (a)  and  (c)  are  given  in  Table  6.1.  In  all  the  cases,  the 
data  from  gauges  at  2.0m  or  4.0m  (remember  these  are  2  experiments  combined) 
were  used  to  synchronize  the  data  with  the  models. 

Figures  6.2  and  6.3  show  comparisons  with  data  from  the  Delft  experiments 
for  case  (a)  of  the  models  WKGS  and  FN4.  Notice  that  at  the  station  5.7m  there 
is  a  phase  mismatch  in  the  data.  This  systematic  error  appears  in  all  the  cases  for 
this  gauge.  Also,  no  data  was  available  for  the  station  23.0m.  Both  the  FN4  and 
WKGS  models  perform  quite  well  for  all  the  gauges  up  to  the  crest  of  the  bar, 
but  as  the  waves  pass  the  back  slope  of  the  bar,  the  WKGS  model  shows  some 
discrepancies  with  the  data.  This  is  due  to  the  aforementioned  decoupling  of  the 
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Figure  6.4:  Linear  dispersion  relationship  as  (nondimensional)  wave  speed  vs. 

wave  frequency.  Present  Model  (dash),  WKGS  (dash-dot),  exact 
(solid).  Dotted  vertical  lines  are  waves  with  periods  Tn  =  (2.02 /n)s. 

higher  harmonics  from  the  primary  longer  wave  which  are  released  as  free  waves 
propagating  with  a  larger  value  of  /i  which  are  more  susceptible  to  inaccuracies. 
The  FN4  model  remains  quite  accurate  even  for  the  gauges  located  after  the 
bar.  To  illustrate  the  inaccuracies  due  to  higher  harmonic  decoupling,  Figure  6.4 
shows  an  alternative  representation  of  the  linear  dispersion  relationship  where  the 
nondimensional  wave  speed  is  plotted  against  the  wave  frequency.  The  vertical 
dotted  lines  indicate  the  location  of  the  frequency  of  the  fundamental  wave  in  case 
(a),  of  which  the  period  is  Ti  =  2.02s,  and  its  harmonics  with  periods  T2  =  7\/2, 
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Ts  =  Ti/ 3,  etc.  Notice  that  the  phase  speed  error  in  the  primary  wave  (Ti)  is  small 
for  both  the  FN4  and  WKGS  models.  As  the  bound  waves  are  released  as  free 
waves  they  travel  with  their  own  speed,  which,  in  the  linear  limit,  are  represented 
by  the  intersection  of  the  vertical  lines  T2,  T3,  etc.  with  each  model’s  dispersion 
curve.  Notice  that  the  errors  in  the  speed  of  the  released  higher  harmonics  starting 
from  T3  for  WKGS  are  considerably  larger  than  for  the  FN4  model. 

Similarly  to  Figures  6.2  and  6.3,  Figure  6.5  shows  plots  for  the  WN4  model 
(this  is  the  present  model,  but  with  the  assumption  0(5)  =  0(/i2)  and  neglecting 
terms  of  0(82fi2, 8(i4, ...).  Apart  from  slight  phase  differences,  the  comparison  is 
about  as  good  as  the  FN4  model,  which  indicates  that  for  this  case,  the  improve¬ 
ment  in  the  dispersion  effects  of  the  FN4  and  WN4  models  over  the  WKGS  model 
is  more  important  than  the  fully  nonlinear  effects  accounted  for  in  WKGS  and 
FN4,  but  not  in  WN4. 

Figures  6.6,  6.7  and  6.8  are  analogous  to  Figures  6.2,  6.3  and  6.5,  but  for 
case  (c)  (see  Table  6.1).  Notice  that  in  this  case  the  incident  wave  has  twice  the 
amplitude  and  about  2/5  of  the  wavelength  of  case  (a).  Before  the  waves  reach 
the  back  slope  of  the  bar,  FN4  and  WKGS  perform  quite  similarly,  although  some 
phase  differences  are  apparent.  Model  WN4  does  not  perform  as  well  in  this 
case  due  to  its  weak  nonlinearity  assumption.  As  the  waves  pass  over  the  bar,  the 
higher  harmonic  decomposition  combined  with  nonlinear  effects  are  strong  enough 
in  this  case  to  make  the  three  models  give  very  different  results,  with  FN4  being 
the  most  accurate,  giving  very  good  agreement  except  for  slight  phase  differences. 
Notice  that  for  this  case,  WN4  does  not  perform  nearly  as  well  as  FN4,  and  also 
qualitatively  worse  than  the  WKGS  (accurate  to  0(fi2),  but  fully  nonlinear)  for 
all  gauges  up  15.7m.  This  result  is  a  strong  indication  that  the  improvement  in 
linear  dispersion  is  not  always  more  important  than  the  fully  nonlinear  effects, 
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Figure  6.7:  Comparisons  of  free  surface  displacement  with  case  (c)  of  Delft  ex¬ 
perimental  data  at  several  gauge  locations.  FN4  (dash-dot),  data 
(solid). 
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contrary  to  the  generalizing  conclusion  of  Dingemans  (1994).  Referring  back  to 
Figure  6.4,  again,  after  the  waves  pass  the  bar,  the  higher  harmonics  are  released 
as  free  waves.  For  case  (c)  the  primary  (incident)  wave  is  indicated  by  the  vertical 
dotted  line  labeled  T2 ,  and  its  second  and  third  harmonics  are  represented  by  the 
even  indexes,  that  is,  T4  and  Tq  respectively.  Notice  that  the  error  in  the  speed 
of  the  primary  wave  is  negligible  for  the  WKGS  model.  In  the  second  harmonic 
(T4)  the  error  for  WKGS  is  fairly  high,  and  for  FN4,  although  not  negligible,  is 
considerably  smaller,  and  the  same  being  the  case  of  the  released  third  harmonic 

(T«). 

Figure  6.9  shows  comparisons  of  the  absolute  value  of  the  amplitudes  of 
the  Fourier  transform  of  one  wave  period  of  the  time  series,  between  both  FN4 
and  WKGS,  and  the  data  points  at  each  gauge  location  for  both  cases  (a)  and 
(c).  Figure  6.10  shows  similar  plots  for  FN4  and  WN4,  where  FN4  results  are 
identical  to  those  in  Figure  6.9.  Also  shown  are  snapshots  of  the  free  surface 
elevation  and  the  position  of  the  bar  (out  of  scale).  In  both  cases  (a)  and  (c), 
the  WKGS  model  tends  to  overpredict  the  higher  harmonics  after  the  crest  of 
the  bar.  For  case  (a)  the  FN4  and  WN4  models  give  very  similar  results,  with 
some  slight  underpredictions  by  WN4  of  the  amplitudes  of  the  released  third 
and  fourth  harmonics  after  the  bar  crest.  In  case  (c)  WN4’s  inability  to  generate 
higher  harmonics  accurately  due  to  the  weak  nonlinearity  assumption  is  evident  in 
the  underprediction  of  the  decomposed  higher  harmonics.  Notice  the  modulation 
present  in  the  fundamental  wave  before  the  bar,  shown  by  all  three  models,  caused 
by  partial  wave  reflection  from  the  front  of  the  bar.  Notice  also  that  for  case  (c)  the 
FN4  model  slightly  overpredicts  the  third  and  fourth  modes  around  the  toe  of  the 
front  face  of  the  bar.  This  is  due  to  numerical  error  introduced  by  the  high  order 
derivative  terms,  which  are  undefined  functions  at  that  location.  When  necessary, 
the  solution  was  filtered  (see  Chapter  4)  to  avoid  high  frequency  contamination 
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Figure  6.9:  Comparisons  of  the  spatial  variation  of  the  Fourier  components  of  the 
free  surface  displacement  with  cases  (a)  and  (c)  of  Delft  experimental 
data.  Bottom  panels  show  the  free  surface  elevation.  FN4  (solid), 
WKGS  (dash-dot),  data  (circles). 
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Figure  6.10:  Comparisons  of  the  spatial  variation  of  the  Fourier  components 
of  the  free  surface  displacement  with  cases  (a)  and  (c)  of  Delft 
experimental  data.  Bottom  panels  show  the  free  surface  elevation. 
FN4  (solid),  WN4  (dash-dot),  data  (circles). 
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Table  6.2:  Index  of  agreement  d,-. 


Gauge 

Location  (m) 

case  (a) 

case  (c) 

WKGS 

FN4 

WN4 

WKGS 

FN4 

WN4 

2.0 

0.998 

0.998 

0.998 

0.997 

0.996 

0.998 

4.0 

0.996 

0.996 

0.996 

0.997 

0.997 

0.984 

10.5 

0.995 

0.995 

0.995 

0.982 

0.986 

0.997 

12.5 

0.999 

0.999 

0.998 

0.997 

0.995 

0.927 

13.5 

0.996 

0.995 

0.987 

0.996 

0.996 

0.990 

14.5 

0.995 

0.997 

0.993 

0.979 

0.971 

0.883 

15.7 

0.995 

0.996 

0.980 

0.973 

0.993 

0.977 

17.3 

0.975 

0.995 

0.972 

0.880 

0.973 

0.934 

19.0 

0.973 

0.982 

0.943 

0.968 

0.987 

0.970 

21.0 

0.927 

0.993 

0.962 

0.948 

0.965 

0.931 

problems.  In  general,  as  in  the  case  of  the  time  series  plots,  the  FN4  agrees  with 
the  data  much  better  than  WKGS  and  than  WN4  for  case  (c). 

Table  6.2  shows  the  index  of  agreement  dt-  defined  by  (6.1)  of  the  models 
FN4,  WN4,  and  WKGS,  with  both  cases  (a)  and  (c)  of  the  Delft  experiments  for 
all  gauges  except  5.7 m  and  23.0m.  Of  course,  the  differences  in  d;  between  the 
models  should  only  have  significance  when  they  are  larger  than  d,-  for  the  incident 
wave  (gauges  2.0m  and  4.0m).  The  results  confirm  that  the  best  performance 
is  from  the  FN4  model,  with  only  one  case  where  WKGS  gave  a  slightly  better 
result  (case  (c),  gauge  14.5m)  due  to  a  slightly  larger  phase  mismatch  in  FN4.  It 
is  clear  that  the  WKGS  model  outperforms  the  WN4  model  around  the  bar  crest 
(gauges  12.5m  through  14.5m),  but  as  the  waves  reach  deeper  water  (importance 
of  nonlinearity  and  dispersion  switch),  WKGS  loses  accuracy.  Although  WN4  has 
much  more  accurate  dispersion  relationship  in  deeper  water  than  WKGS,  since 
it  was  not  capable  of  generating  higher  harmonics  properly  while  the  waves  were 
shoaling,  the  overall  solution  becomes  inaccurate  after  the  bar.  This  confirms  the 
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Figure  6.11:  Sketch  of  wave  flume  of  the  Ohyama  experiment.  All  dimensions 
in  (m) 


importance  of  the  full-nonlinearity  assumption  made  in  the  WKGS  derivation  but 
not  in  the  WN4. 


6.2  The  Ohyama  Experiment 

In  this  section  we  show  comparisons  of  the  models  FN4  and  WKGS  model 
with  the  experiment  by  Ohyama  et  al.  (1994)  (referred  here  as  simply  the  Ohyama 
experiment).  Computations  with  the  WN4  model  were  not  performed  for  this  case. 
A  sketch  of  the  wave  flume  with  the  gauge  locations  is  shown  in  Figure  6.11.  We 
now  summarize  the  experimental  setup.  The  wave  flume  is  65 m  long  and  1.0m 
wide.  The  total  depth  of  the  flume  is  1.6m.  The  location  of  the  center  of  the  bar 
was  28.3 m  from  the  piston-type  wavemaker.  All  other  relevant  dimensions  can 
be  seen  in  Figure  6.11.  The  measurements  were  performed  before  the  point  when 
waves  reflected  from  the  bar  reached  the  wavemaker.  At  the  right  end  of  the  flume, 
waves  were  absorbed  by  the  presence  of  coarse  materials  to  dissipate  the  energy. 
A  total  of  six  tests  were  performed  with  three  different  incident  wave  periods 
(1.34s, 2. 01s, 2. 68s)  each  for  two  different  wave  amplitudes  (0.0125m, 0.025m).  No 
wave  breaking  occurred  in  any  of  the  tests.  The  data  was  obtained  by  digitization 
of  the  plots  (performed  by  Andrew  Kennedy)  from  the  original  article.  The  only 
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Table  6.3:  Incident  wave  characteristics  for  the  Ohyama  experiment 


Case  (2) 

Case  (4) 

Case  (6) 

Wave  amplitude  (m) 

0.025 

0.025 

0.025 

Wave  period  (s) 

1.341 

2.012 

2.683 

-e 

III 

1.299 

0.769 

0.555 

5  =  ao//io 

0.050 

0.050 

0.050 

time  series  available  for  comparisons  were  the  ones  at  station  3  and  5,  for  all  three 
wave  periods,  and  the  highest  of  the  two  amplitudes  (0.025m).  Fourier  amplitudes 
were  available  for  the  same  wave  conditions  but  at  all  measurement  stations.  Time 
series  were  synchronized  at  station  3.  We  refer  to  the  three  tests  as  cases  (2),  (4), 
and  (6),  as  in  Ohyama  et  al.  (1994).  The  incident  wave  conditions  are  summarized 
in  Table  6.3.  The  incident  wave  conditions  are  similar  in  the  Ohyama  and  Delft 
experiments.  The  major  difference  between  the  two  experiments  is  that  the  bar 
in  the  Ohyama  experiment  is  much  shorter  and  with  much  steeper  slopes  than 
the  one  in  the  Delft  experiments,  more  reminiscent  of  a  submerged  rubble  mound 
structure.  The  steep  slopes  add  extra  difficulty  for  the  models’  performance,  since: 
(i)  the  models’  dispersion  properties  are  optimized  assuming  constant  depth;  (ii) 
the  assumption  that  the  vertical  velocity  is  0(fi2)  times  the  horizontal  velocity 
is  violated  at  steep  slopes.  Smoothing  of  the  corners  of  the  bar,  besides  filtering 
every  100  time  steps  was  necessary  to  prevent  spurious  high  frequency  noise  to 
contaminate  the  solutions.  To  smooth  the  corners  of  the  bar  we  applied  a  3-point 
average  by  Shapiro  (1970)  five  times.  Since  the  waves  are  progressively  longer 
from  case  (2)  through  (6),  we  expect  that  the  Boussinesq  models  will  perform 
best  in  case  (6),  and  worst  in  case  (2).  We  also  expect  higher  mismatches  between 
models  and  data  at  station  5  than  at  station  3,  due  to  increasing  errors  in  the 
phase  of  the  decomposed  higher  frequency  bound  waves  as  they  reach  the  deeper 
water  behind  the  bar. 
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Figures  6.12,  6.13,  and  6.14  show  comparisons  of  the  FN4  and  WKGS 
models  with  data  for  cases  (2),  (4),  and  (6),  respectively.  Notice  that  for  cases  (2) 
and  (4)  the  FN4  model  shows  a  mismatch  in  the  phase  speed  at  station  5,  and  an 
underprediction  of  the  wave  crests  and  troughs,  an  indication  that  even  the  fully 
nonlinear,  0(ji4)  model  has  limited  ability  to  predict  waves  past  a  submerged  bar 
with  very  steep  slopes,  if  the  waves  are  short  enough.  For  case  (6)  the  FN4  model 
agrees  very  well  with  the  data.  For  all  three  cases,  the  WKGS  model  has  poor 
qualitative  agreement  with  the  data  at  station  5,  mostly  due  to  phase  errors  and 
overprediction  of  higher  harmonics  behind  the  bar. 

Figures  6.15,  6.16,  and  6.17  show  comparisons  of  the  Fourier  amplitudes 
along  the  flume  between  both  FN4  and  WKGS,  and  the  data  points  at  each  station 
for  cases  (2),  (4),  and  (6)  respectively.  In  all  cases,  the  models  predict  well  the 
Fourier  amplitudes  before  the  back  face  of  the  bar.  For  case  (2),  WKGS  gives 
slightly  better  prediction  of  the  second  harmonic  at  stations  4  and  5  than  FN4, 
but  once  again  strongly  overpredicts  the  third  and  fourth  harmonics  at  those 
stations.  For  case  (4),  the  FN4  model  gives  better  prediction  than  WKGS  for 
all  but  the  third  harmonic,  which  WKGS  agrees  slightly  better  with  the  data. 
For  case  (6),  both  models  agree  reasonably  well  with  the  data,  with  FN4  having 
a  better  prediction  of  the  third  harmonic  at  station  5  and  the  WKGS  model 
matching  the  fourth  fourth  harmonic  slightly  better  at  that  same  station.  For 
this  case,  the  deviations  from  the  data  in  the  time  series  computed  by  WKGS  at 
station  5  are  probably  due  to  phase  errors,  which  is  not  detected  by  the  Fourier 
amplitudes  comparisons. 

Similarly  to  the  case  of  the  Delft  experiment,  table  6.4  shows  the  index  of 
agreement  between  the  models  WKGS  and  FN4,  and  the  data  from  the  Ohyama 
experiment  for  cases  (2),  (4),  and  (6),  stations  3  and  5.  Notice  that  for  cases  (2) 


101 


(iu)h.  (ui)  U 


Station  3  -  (a)  Station  5  -  (b) 


Figure  6.12:  Comparisons  of  free  surface  displacement  with  case  (2)  of  the 
Ohyama  experimental  data  at  stations  3  and  5.  FN4  (upper  panels 
-  a,b),  WKGS  (lower  panels  c,d),  data  (circles). 


Table  6.4:  Index  of  agreement  dt. 


Station 

case  (2) 

case  (4) 

case  (6) 

WKGS 

FN4 

WKGS 

FN4 

WKGS 

FN4 

3 

0.994 

0.998 

0.991 

0.991 

5 

0.921 

0.914 

0.927 

mi 

0.945 

0.976 
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Figure  6.13:  Comparisons  of  free  surface  displacement  with  case  (4)  of  the 
Ohyama  experimental  data  at  stations  3  and  5.  FN4  (upper  panels 
-  a,b),  WKGS  (lower  panels  c,d),  data  (circles). 


i\(m)  r \(m) 


Station  3  -  (a) 


Station  3  -  (c) 


Station  5  -  (b) 


Figure  6.14:  Comparisons  of  free  surface  displacement  with  case  (6)  of  the 
Ohyama  experimental  data  at  stations  3  and  5.  FN4  (upper  panels 
-  a,b),  WKGS  (lower  panels  c,d),  data  (circles). 
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(w)  1  U|  (ui)  fill  (ui)fh\  (Ui)\\  (Ui) u 


Fourier  amplitudes,  case  (2) 


x  10~2 


Figure  6.15:  Comparisons  of  the  spatial  variation  of  the  Fourier  components  of 
the  free  surface  displacement  with  case  (2)  of  the  Ohyama  experi¬ 
mental  data.  Bottom  panel  shows  the  free  surface  elevation.  FN4 
(solid),  WKGS  (dash-dot),  data  (circles). 
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T\(m)  hi  I  (m)  lr|  I  (mj  hi  I  (m)  hi {\(m) 


Fourier  amplitudes,  case  (4) 


x  10'2 


Figure  6.16:  Comparisons  of  the  spatial  variation  of  the  Fourier  components  of 
the  free  surface  displacement  with  case  (4)  of  the  Ohyama  experi¬ 
mental  data.  Bottom  panel  shows  the  free  surface  elevation.  FN4 
(solid),  WKGS  (dash-dot),  data  (circles). 
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x  jq-2  Fourier  amplitudes,  case  (6) 


Figure  6.17:  Comparisons  of  the  spatial  variation  of  the  Fourier  components  of 
the  free  surface  displacement  with  case  (6)  of  the  Ohyama  experi¬ 
mental  data.  Bottom  panel  shows  the  free  surface  elevation.  FN4 
(solid),  WKGS  (dash-dot),  data  (circles). 
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and  (4),  the  results  indicate  a  better  agreement  with  the  data  by  WKGS  than  by 
FN4.  By  inspecting  time  series  comparisons  in  Figures  6.12  and  6.13,  it  is  clear 
that  the  better  agreement  index  for  the  WKGS  model  is  only  due  to  a  systematic 
phase  error  by  the  FN4  model,  which,  overall  has  a  better  qualitative  agreement. 
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Chapter  7 


CONCLUSIONS  AND  RECOMMENDATIONS 


7.1  Conclusions 

A  Boussinesq-type  model  with  0(1)  nonlinearity  and  0(/j,4)  dispersion  and 
vertical  dependence  was  developed.  By  conveniently  defining  one  of  the  dependent 
variables  as  the  weighted  average  of  the  velocity  potential  at  two  distinct  water 
depths,  it  is  possible  to  achieve  an  accurate  (4,4)  Pade  approximant  form  for  the 
linear  dispersion  relationship.  A  major  improvement  over  the  existing  second  order 
models  has  been  found  in  the  prediction  of  the  linear  internal  flow  kinematics. 

A  perturbation  approach  was  carried  out  to  analyze  random  wave  second 
order  nonlinear  interactions  and  it  has  been  shown  that  the  FN4  predicts  very  well 
the  transfer  coefficients  of  super  and  subharmonics  generation  over  a  wide  range 
of  water  depths.  The  model’s  cubic  nonlinear  Schrodinger  equation  governing  the 
propagation  of  wave  group  envelope  was  obtained  by  a  standard  WKB  perturba¬ 
tion  multiple  scales  approach  and  its  coefficients  were  compared  to  those  of  the 
full  model  as  well  as  of  WKGS  0(//2)  model.  The  model’s  cubic  term  coefficient 
was  shown  to  have  good  agreement  with  the  full  model. 

A  numerical  implementation  of  the  1-dimensional  version  of  the  model  was 
used  to  simulate  wave  evolution  over  arbitrary  bottom  topography.  The  numerical 
model  included  absorbing  sponge  layers  to  simulate  radiation  boundary  conditions, 
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and  generation  of  waves  inside  the  domain  by  the  inclusion  of  a  source  function  in 
the  system  of  equations.  The  source  function  was  derived  from  Green’s  functions 
linear  theory,  but  it  was  shown  that  it  also  works  well  for  finite,  relatively  small 
amplitude  waves. 

Several  numerical  computations  were  carried  out  for  solitary  waves  propa¬ 
gating  over  both  horizontal  and  sloping  bottoms,  and  many  solitary  wave  proper¬ 
ties  were  compared  with  both  the  exact  solution  and  other  models.  It  was  shown 
that  the  FN4  model  agrees  better  with  the  exact  solution  than  does  the  WKGS 
model.  Comparisons  between  FN4  and  GN3  (third  order  Green-Naghdi  model) 
with  the  exact  solution  showed  that  GN3  agrees  better  than  FN4  for  some  prop¬ 
erties.  For  sloping  bottom,  FN4  and  WKGS  models  were  compared  with  the  very 
accurate  BEM,  and  FN4  agreed  with  BEM  better  than  WKGS.  A  set  of  GN3 
evolution  equations  for  variable  depth  is  not  available  at  this  time,  therefore  no 
comparisons  between  FN4  and  GN3  are  available  for  this  case. 

Finally,  numerical  computations  of  FN4  and  WKGS  were  compared  to 
several  laboratory  measurements  of  waves  propagating  over  submerged  bars,  and 
FN4  generally  gave  better  agreement  with  the  data.  The  weakly  nonlinear  version 
of  the  present  model,  WN4,  was  also  compared  to  some  of  the  experimental  data. 
The  results  showed  that  the  higher  order  nonlinear  terms  neglected  in  WN4  are 
essential  for  accurate  prediction  of  higher  amplitude  waves,  specially  if  these  waves 
are  in  deeper  water,  since  the  higher  harmonics  are  more  susceptible  to  phase 
errors. 

7.2  Recommendations  for  Future  Work 

The  numerical  implementation  of  the  present  model  was  found  to  have  the 
following  minor  problems:  (i)  for  very  steep  (nonlinear)  waves,  the  convergence 
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of  the  iterative  time  integration  can  be  fairly  slow,  besides  that,  in  some  cases,  a 
very  low  under-relaxation  parameter  is  needed  to  avoid  divergence  of  the  solution. 
Further  investigation  is  required  to  improve  the  model  in  this  aspect,  (ii)  A  better 
treatment  of  the  derivatives  of  the  water  depth  at  locations  where  they  are  singular 
could  improve  the  performance  of  the  model  in  that  as  it  is  now,  filtering  has  to 
be  applied  to  the  solution  every  several  time  steps,  and  in  cases  of  very  abrupt 
transition,  the  sharp  corners  of  the  bottom  needs  to  be  smoothed. 

In  this  work  we  used  Green’s  function  theory  (linear)  to  derive  the  necessary 
source  function  for  wave  generation  inside  the  domain.  In  cases  of  relatively  high 
amplitude  waves,  the  waves  generated  by  the  source  function  are  not  a  solution  to 
the  equations,  and  therefore,  undesirable  frequencies  may  appear  in  the  solution 
near  the  source  region  to  “make  up”  for  the  differences.  Further  investigation  is 
needed  in  this  direction  to:  (i)  quantify  the  validity  of  the  use  of  a  source  function 
derived  from  linear  theory  with  a  nonlinear  set  of  equations,  (ii)  derive  a  source 
function  which  generates  solutions  of  the  nonlinear  problem  for  regular  waves  such 
as  Stokes  and  cnoidal  waves. 

Lastly,  the  numerical  implementation  can  be  expanded  to  two  spatial  di¬ 
mensions  in  a  straightforward  manner,  although  the  computational  cost  involved 
in  solving  the  system  would  be  quite  high  for  present  standards.  This  step  is  also 
left  for  future  work. 
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Appendix  A 


APPLICATION  OF  NWOGU’S  METHOD  AT  O(^) 


The  procedure  of  Nwogu  (1993)  rests  on  choosing  the  potential  or  velocity 
at  an  elevation  za  in  the  water  column  such  that  the  resulting  linear  dispersion 
relationship  of  the  model  is  optimal  in  some  sense.  The  dispersion  relationship  in 
a  model  retaining  terms  to  0(fi2)  is  given  by 


1  -  (a  +  l/3)/U2 

1  —  ajj, 2 


(A.l) 


where 

«  =  \zl  +  (A-2) 

The  choice  a  =  —2/5  reproduces  the  (2,2)  Pade  approximant,  while  the  choice 
a  =  —0.39  proposed  by  Nwogu  minimizes  the  error  (in  a  least-square  sense)  in 
the  dispersion  relation  over  the  range  0  <  fx  <  n. 


Following  this  procedure,  we  extend  the  approximate  expression  for  the 
velocity  potential  (in  terms  of  </>a)  to  0(fi4)  and  obtain 

0  —  0a  +  [(1  +  £«)2  “  (1  +  ^)2]  V20<* 

+  ~  [5(1  +  zQ)4  —  6(1  +  ^a)2(l  +  z)2  +  (1  +  ^)4]  V2V20a 
+  0(f)  (A.3) 
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This  expression  is  used  in  linearized  versions  of  (2.4)  and  (2.6)  to  obtain  the  linear 
model 

Vt  +  VVa  +  ^2(a  +  l/3)V2VVa  +  /(5a2+4a  +  4/5)V2V2VVa  =  0 

rj  +  <j>at  +  yu2aV2</>at  +  /^(2  +  5a)V2V2</>ai  =  0  (A.4) 

0 

The  corresponding  linearized  dispersion  relation  is  given  by 


ui2  = 


1  —  (a  +  1/3  )fj,2  +  (5a2  +  4a  +  4/5  )fi4 


(A.5) 


1  —  an2  +  (5  a2/6  +  a/3  )/j,4 

This  is  equivalent  to  Nwogu’s  result  if  terms  of  0(p4)  are  dropped.  The  resulting 
dispersion  relation  contains  only  a  single  parameter,  and  there  is  no  choice  of  a 
which  reproduces  the  desired  (4,4)  Pade  approximant: 


a;2  =  1  +  (V9)^2  ±  (1/945 )jS_  + 


1  +  (4/9)yu2  +  (1/63  )fi4 


(A-6) 


It  has  been  shown,  therefore,  that  by  using  Nwogu’s  method  of  choosing 
the  potential  at  a  depth  zQ  as  the  dependent  variable  in  a  4th  degree  polynomial 
in  z,  it  is  not  possible  to  obtain  the  desired  (4,4)  Pade  approximant. 
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Appendix  B 


DERIVATION  OF  THE  SCHRODINGER  EQUATION 
ASSOCIATED  WITH  THE  PRESENT  MODEL 

Here  we  present  the  derivation  of  the  Schrodinger  equation  for  the  present 
model  in  dimensional  form.  The  nondimensional  final  results  are  presented  in 
Chapter  3. 

The  constant  depth  versions  of  the  evolution  equations  for  4>  and  77  in 
dimensional  form  are: 

m  +  v-{ff  [v<A+|(B-|ff2)w2<A 

+  i  (b2-\d-  ±BH2  +  iff4)  VV2vV] }  =  0  (B.l) 

m  +  k  +  \  |v^f  +  \  (b  -  ff2)  V*i, 

+  -  (b2  —  -D  -  BH 2  +  iff4!  V2V24>, 

4  V  6  6  / 

+  v<£-v[i(s-tf2)  v2<£ 

+  ^  (b2  -  -  BH 2  +  i H 4)  V2V2<£] 

+  iv|(B-F2)2(v2^)|2  +  ^2(v2^)2 

+  ^  (bH2  -  V2</>V2V2</>  =  0  (B.2) 
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Assuming  8  as  our  small  parameter,  we  expand  for  our  dependent  variables  as: 


V  =  <%i  +  8\  +  83r)3 

<j> ==  8(f>  i  +  82  <f>2  +  53<^3  (B.3) 

and  introduce  multiple  scales  for  the  independent  variables  as: 

t  =  t'  +  8?  +  8H'  =  t'  +  T4+T2 

x  =  xf  8x‘ '  4“  82xf  =  x1  -f-  Af  4~  -A2 

V  =  8y'  +  82y'  =  Y1+Y2  (B.4) 

Before  we  proceed,  we  define  the  following  operators: 

-bl(')  =  (')tt  4“  S,h(')xx  fjCl h  (')xxxx  4“  yC2h  (')xxxxxx 

4”  (')xxtt  C 4ll  (* fxxxxxx  (B*5) 

L2(-)  =  ~2(-)tTl  +  2gh(-)xx1  -h  2gh(-)xX!  +  4C1h3(-)xxxX1 
4-  6gC2h5 (-jxxxxxX}  4-  2C3h2(-')XxtTi 

+  2 C3h2(-)xXltt  ~  2C4h4(-)xxxxtTl  -  4 C4h\-)xxxXltt  (B.6) 

Lz{-)  =  -(OtiTi  -  2(-)iT2  4-  gh(-)x \Xi  4-  2gh(-)xx2  4-  gh{-)y ,y, 

—  C\h  ^  4“  4( * )xxxX2  4-  2(-)xa;y1  y, ] 

4"  C2h  \Yi('')xxxxXiXi  4"  Q(’)xxxxxX2  4"  3( •  )xia;xy]  yj  ] 

4-  C3h5  [(-)xxTiTi  4-  2(-)xxtT2  +  4(-) xXitT 1 
4-  (-fxiXitt  4-  2(-)xx2tt  +  (OviVitt] 

C4h  [(')xxxTi  T\  4"  ( ’ )xxxxtT2  4"  2(-^XXXxtT2  4"  8(')xxxX]  tT\ 

4-  6(-)a;a:X1X1tt  4-  4(-)xxxx2tt  4-  2(-)xa;y1yltt]  (B.7) 

^i(’)  =  (’)t  —  C3h2(-)xxt  4-  C4h4{-)xxxxt  (B.8) 

L2{-)  =  (-)t!  -  C3h2  [(OxxTj  4-  2(-)xXlt] 
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(B.9) 


+  ^4^  \^*)xxxxT\  “h  ('JacafttAif] 

L'A')  —  {-)t2  ~  C3h2  [2(-)xXiTi  +  {‘)xxT2  +  2(‘)xX3t  +  (■)x1x1t  +  (Onvit] 

+  C4hA  [{•)XxxxT2  +  4(-)mXiTi  +  ^(■)xxxX2t 

+  Gi^xxX^!  t  +  2(-)x;ry1y1<]  (B.IO) 

where  C i,  C2,  C3,  C4  are  given  in  the  section  below.  For  each  order  n  =  0, 1, ...  we 
assume  the  following  solution  corresponding  to  a  nonlinear  wave  train  propagating 
mainly  in  the  x  direction,  but  allowing  the  amplitudes  (f)nm  and  r\nm  to  have  slow 
variations  in  space  and  time: 

Vn  =  £  rlnm(X1,X2,Y1,Y2,T1,T2)Em  (B.ll) 

m——n 

<f>u  =  £  cf>nm(X1,X2,YuY2,TuT2)Em  (B.12) 

m——n 

where  E  =  e^kx'  -wt\  n  indicates  the  order  of  the  solution,  and  negative  m  indices 
indicate  the  complex  conjugate  of  —m  for  both  gnm  and  <f>nm.  Substituting  (B.3) 
and  (B.4)  into  the  t  derivative  of  (B.2)  minus  g  times  (B.l),  and  into  (B.2),  and 
ordering  the  problem  in  powers  of  S,  we  obtain  for  n  =  1: 

Lx<t>x  =  0  (B.13) 

ggi  +  L\<$>  1  =  0  (B.14) 

Substituting  (B.ll)  and  (B.12)  with  n  =  0  into  (B.13)  and  (B.14),  we  obtain: 

DX(f>n  =  0  (B.15) 


and 


<t>  11  = 


-igA 


2w[i  +  c3(khy  +  c4(khyy 


A  =  2?7ii 


(B.16) 


116 


where 


Dx  =  —gk(kh)  [l  +  C^khf  +  C2{kh)4]  +  u2  [l  +  C3(kh)2  +  C4(kh)4]  (B.17) 

and  (B.15)  is  the  linear  dispersion  relationship.  We  arbitrarily  assumed  that 
7710  =  0.  We  now  proceed  to  seek  a  “slowly  varying  equation”  for  A. 

For  n  —  2,  we  obtain  the  equations: 

Li<f>2  +  L2(f>i  +  g{rfr(f>ix)x  —  gC3h2(r]14)ixxx)x 

+  gCAh4{q4(j>Xxxxxx)x  +  h(g1(f>lxxt)t  -  <f>ix4>ix)t 

—  Cih3(gi(f)ixxxxt)t  +  C3hz{(j)ix^  ixxx)t 

C4h  {^<j)ix(j)\ xxxxx)t  2^3^  (4*1  xxx&l xxx)t 

~  2^2^lxx^>lxx^t  4(4)lxx(t)lxxxx)t  —  0  (B.18) 

and 

912  +  L[(j)2  +  L'2(f>i  —  hr)i4>ixxt  +  ~<j>\x4>\x  +  Cih3gi<f>ixxxxt 

+  C3h2  4>ix(f>ixxx  +  C4h44>ix4>ixxxxx 

+  2^2<^1xx<^1xx  -  Cih24>ixx4>ixxxx  =  0  (B.19) 

Substituting  (B.ll)  and  (B. 12)  with  n  =  1  into  (B.18),  and  matching  the  coeffi¬ 
cients  of  each  power  of  E ,  we  obtain  (all  terms  ~  E°  are  zero  for  this  equation  at 
this  order): 

(i)  terms  ~  E1: 

2ico  [l  +  C3{kh)2  +  C4(kh)4]  [<f>UTi  +  Cg(f>nXl]  =  0  (B.20) 
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where  Cg  =  dto/dk,  from  the  dispersion  relationship  (B.15). 


(ii)  terms  ~  E 2: 

D2<f>22  -  2gk2  {l  +  C3{kh)2  +  C4{kh)4 

-  y  [l  +  CxWjjwn 

-  iuk 2  [l  +  2 C3(kh)2  +  4C4(kh)4 

+  Cl(kh)4-(kh)2-2C1(kh)4]<f>211  (B.22) 

where 

D2  =  -4gk{kh)  [l  +  4Cx{kh)2 +  %C2(kh)4] 

+  4a;2  [l  +  4 C3{kh)2  +  8 C4(kh)4]  (B.23) 

Substituting  (B.16)  into  (B.22),  we  obtain  after  some  algebra: 

fa  =  ^^-^[l  +  c^khy  +  c.ikh)4}2 

-  2(khf  [l  +  Cxikh)2  +  C2(kh)4]  [l  +  Cl{kh)2] 

+  1  +  (2 C3  -  1  ){kh)2  -  2 C4(kh)4  +  Ci(kh )4  -  Cxikh)4}  {( kh )2 
+  15(C4  -  C2)(kh)4  +  12(C'iC4  -  C2C3)(kh)6'^~1  (B.24) 

Substituting  (B.ll)  and  (B.12)  with  n  =  1  into  (B.19),  and  matching  the  coeffi¬ 
cients  of  each  power  of  E ,  we  obtain: 

(i)  terms  ~  E° : 

gg 20  =  ^iotj  ~  [(2C3  +  l)(kh)2  +  l|  A;2<^>h<5^i_i 

+  —  rjnfii-i)  ~  (2Ci  +  C\  +  2C±){kh)Ak2(j)ii(f)i-i 

+  icoCiK4 h3 -  7711^1— i )  (B.25) 
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Using  (B.16)  into  (B.25),  we  obtain  after  some  algebra: 


9  20  —  ^ioj, 
9 


+ 


-  [(2 C3  +  l){khf  +  1]  +  [—2 C3  -  C2  -  2 C4\  ( khy 
1  +  Cx{khf  +  C2{kh)4 
|A|2 


+  2{khf  +  2C1{kh)4} 


Ah  [1  +  C3(kh f  +  C5(kh)4] 


(B.26) 


(ii)  terms  E1: 

gi i2i  =  [  1  +  C3(kh)2  +  C4{kh )4]  hi  -  [l  +  C3{khf  +  C4{kh)4} 

+  2 ukh2  [C3  +  2C4(kh)2]  <f>  11Xl  (B.27) 


Since  <j> 21  can  be  absorbed  into  0n,  we  set  <j> 21  =  0  without  loss  of  generality. 
Substituting  (B.16)  into  (B.27),  we  obtain  after  some  algebra: 

C3{khf  +  2  C4(kh)4 


V21 


^—A 

2uj 


T\ 


'k[i  +  c3(khy  +  c4{khy 


(B.28) 


(ii)  terms  ~  E 2: 


9922  —  2zu;  [l  +  AC3(kh )2  +  16C'4(A:/i)4J  <^22 

d-  iujk2 -(-  — fc202i  -j-  iu)Cik4h3r}n(f>ii  -1-  C3k4h2<p 2^ 

+  C4k6h44>2n  +  \clk6h4(f)2n  -  \k4h2<j>\x  -  Cxk3h44> 21  (B.29) 

z  z 


Substituting  (B.16)  and  (B.24)  into  (B.29),  we  rewrite  (B.29)  is  the  schematic 
form: 


V22 


(B.30) 


where  S22  is  a  complicated  function  of  kh. 
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We  are  now  ready  to  move  on  to  the  next  order,  n  =  3.  At  this  order,  a 
great  amount  of  algebra  is  necessary  in  the  evaluation  of  all  the  nonlinear  terms, 
and  since  we  already  have  evaluated  all  the  dependent  variables  appearing  in  these 
terms  from  the  previous  orders,  from  this  point,  we  only  outline  the  steps  before 
we  give  the  final  result.  For  n  =  3,  the  combined  equation 

djdt  (B.2)  -  gx  (B.l)  becomes: 

Lifa  +  L2&2  +  L>3fa  +  %i  +  +  Z3  =  0  (B.31) 

where  Z3  are  the  nonlinear  cubic  terms  involving  products  of  rji  and  fa,  and 
their  fast  time  ( t )  and  spatial  (x)  derivatives.  Z2  are  the  nonlinear  terms  involv¬ 
ing  products  between  either  771  and  fa  or  772  and  fa ,  and  Z\  are  nonlinear  terms 
involving  products  between  771  and  fa  and  containing  either  a  T\  or  an  X\  deriva¬ 
tive.  We  now  substitute  (B.ll)  and  (B.12)  with  n  =  2  into  (B.31).  To  obtain 
the  Schrodinger  equation  for  A,  only  information  ~  E°  and  ~  E1  is  needed.  The 
only  linear  term  in  (B.31)  that  contributes  to  ~  E°  and  ~  E 1  is  L$fa.  Using  the 
relations 


faiTx  ~  CgfalXx 

(B.32) 

falTiTi  =  Cg^llXiXi 

(B.33) 

=  ~CgfaiXlXl 

(B.34) 

after  some  algebra, 

the  linear  term  Lzfa  becomes: 

L3fa  = 

2 ioj  [l  +  C3(kh)2  +  C4{kh)4}  x 

(^falT2  +  CgfalX2  ~  — ^1  —  E1 

+ 

[— faoTiT!  +  9k  [faox1xl  +  ^lo^y,)] 

(B.35) 
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Where  u>"  is  obtained  by  twice  differentiating  u)  with  respect  to  k  in  the  dispersion 
relationship. 


Substituting  the  relations  between  7711,  7720,  722,  0n>  022,  and  A  in  all  ~  E 1 
nonlinear  terms  in  (B.31),  and  combining  them  with  the  ~  E1  terms  in  (B.35), 
we  find  the  following  equation: 


?  1 C 

At2.+  CgAx 2  -  AXlxx  -  +  ia  1 1^4|2  ^4. 


+  iSfi  <  ^>10x1  + 


1  +  C^g2  +  C47U4  —  g2u> 2  (1  +  Cxg2) 

2lo  (1  +  Cs/j,2  +  C4IJ4) 


<hoXi  >  =  0  (B.36) 


where  o\  is  given  in  nondimensional  form  in  the  next  section. 

Similarly,  we  combine  the  linear  and  nonlinear  ~  E°  terms,  and  after  some 
algebra,  we  obtain 

2 

—0io x,Ti  +  9h  (^loXiXi  +  010k, y, )  =  I -A  lx,  +  So\A\^  (B.37) 


where 

_  _w2  {1  _  ci(kh)2}  [1  +  (1  +  2 Cz){kh)2  +  (2C4  +  C2  +  2C1)(kh)4] 

0  4(kh)2  [1  +  Ci(kh)2  +  C2(kh)4]  1  } 


Now  we  introduce  the  changes  of  variables 


t  =  Xx-  CgTU  r  =  5Tl 


(B.39) 


Substituting  (B.39)  into  (B.37)  and  integrating  once  with  respect  to  £  we  get: 

g2k/(2u))  —  CgSo .  ,  l2 


01Of  — 


C2  -  gh 


-|A|2  +  7o(r) 


(B.40) 
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where  70 (r)  is  an  integration  constant  which  vanishes  for  a  wave  train  beginning 
from  rest  where  A  and  </>10£  tend  to  zero  as  £  -»  00  (see  Mei,  1989).  Substituting 
(B.40)  with  (B.39)  into  (B.36)  with  T2/5  replacing  Ti,  we  have: 

Ar  -  -  Y^-Aym  +  ia  \A\2  A  +  71  {j)A  =  0  (B.41) 

where  71  (r)  absorbs  7o(r).  Without  Yi  derivatives  and  without  the  last  term 
containing  71  (t),  equation  above  is  the  cubic  Schrodinger  equation.  The  last 
term  can  be  absorbed  into  a  new  dependent  variable  by  the  transformation  A'  = 
Ae’f71  dr.  The  nondimensional  version  of  equation  above  is  given  in  Chapter  3. 
The  coefficient  a  =  0\  +  <r2  is  given  in  nondimensional  form  in  the  next  section. 

B.l  Expressions  for  Cubic  Term  Coefficient 


<Tl  = 

+ 

+ 

<X2  = 


— —  —  4  +  16Ci/j-2  —  u>  fl  +  4C3/U2  +  I6C4//4) 
uQiQ?  L  v  n 

2Q- (E20  +  E22)  +  Cifi3  —  u>  2/i  1QiJ 

[1  +  (2  +  5 C)  +  (IOC,  +  17C,  +  4 Cl)  /] 

I1  +  c'’'j2  _  0  +  Cl1’2)] 

^Q2  [4  +  (8C3  +  1/6) /i2] 


(&  ~  ?.<*>) 

2  (Cl  -  l) 


ZqSQi~^2  ('  +  <?./)} +  (1  +  ^) 


Where 


Q\  —  1  +  +  C*4/i4 

Q2  =  4  (l  +  u>  2 j  ^1  4-  4Ci/i2  +  4C4/i4) 
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~  2Ql  ^  ~~  CiV2}  +  4u2Q2jj2  I1  +  (2^3  +  !)  ^  +  (2C*  +  Cl  +  2Ci)  /i4} 

P22  =  -2[Qi-^2(l  +  C'1/u2)]-gr1[2Qi-l-^2-(2C1-(7|)^4] 

E20  =  ^(1  +  Cl^V^^[2^-1  +  ^+(2^  +  C32)/] 

^22  =  g j4'2Qi  P22  (X  +  4Cs^2  +  16°4^ 

+  i|r(1+ Ci^2)  “  [2Qi  - 1  -  + (2Ci  -  c32)  v4} 

and 

Ci  =  -\{B-  1/3);  C2  -  i  (B2  -  5/3  -  D/6  +  1/30) 

C3  =  ~\(B- 1);  C4  =  i  (B2  -  B  -  D/6  +  1/6) 

The  corresponding  crj  and  cr2  for  the  full  boundary  value  problem  are  given  by: 

cosh  4/j  +  8  —  2  tanh2  /2 
16sinh4/2 

*  =  . r±  Y(i+  C»,V 

2u;  (Cg  —  l)  V  2a;  cosh2  ji ) 
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Appendix  C 

DERIVATION  OF  SOURCE  FUNCTION 

Here  we  derive  the  x-direction  source  function  for  regular  waves.  The 
linearized  versions  of  the  mass  and  momentum  equations  for  4>  over  a  fiat  bottom, 
including  the  source  function  is  given,  in  dimensional  form,  as: 

r)t  -  W2^  +  C'1h3V2V2^  +  C2h5V2V2V20=/s(x,J/,t),  (C.l) 

4>t  +  gr]  +  hV2]>-C3h3V2V24>t  +  C4h5\72V2V24>t  =  0,  (C.2) 

where  coefficients  C\,  C2 ,  C3,  C4  are  as  defined  in  Appendix  B.  Taking  the  t 
derivative  of  the  momentum  equation  and  eliminating  rj  from  (C.l)  and  (C.2), 
gives: 

4>tt  -  ghV24>  +  Cigh3V2V2]>-C2gh5V2V2V2<f> 

-  C3h3V2V2^t  +  C4h5V2  V2V2^t  =  -gfs(x,  y,  t).  (C.3) 

We  introduce  the  following  transformations: 

fa,y,t)  =  (C-4) 

f(x,y,t)  =  (1-)  }(x)e'Xye-'“,d\du.  (C.5) 

Substituting  (C.4)  and  (C.5)  into  (C.3)  we  have: 

a$G]  +  b^[4]  +  c$2]  +  4  =  gf ,  (C.6) 
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where  the  numbers  in  brackets  denote  order  of  x  derivatives,  and 


a  =  C2gh 5, 

b  =  -Cigh3  +  C4h4u2-3C2gh5X2, 

c  =  gh-C3h2u2 +2Cigh3X2 -2C4h4X2u2 +  ZC2gh*X4, 

d  =  u2  -gh\2  +  C3h2X2u2  -Cigh3X4  +  C4h4X4u2  -C2gh3X6.  (C.7) 


Now  we  multiply  (C.6)  by  a  Green’s  function  G(£,  x),  and  integrate  the  product 
with  respect  to  £,  from  — oo  to  +oo,  which  gives: 


[  °°  (aG[6]  +  bG[4]  +  cG[2]  +  dG)  <£d£ 

+  a  [g/1  -  GW^41  +  Gmf]  -  <*P  +  d4'/1  +  G^ 

Of1  -  G«f]  +  G'2'/'  -  Gl3^l+“ 

-  —  OO 

G^1]  -  G^]+°°  =  g  f+°°  GM, 

J  —00  J  —00 


1  +00 


+  b 

+  c 


(C.8) 


where  the  numbers  in  brackets  denote  order  of  £  derivatives.  Notice  that  £  is 
a  dummy  variable  and  x  is  now  an  arbitrary  fixed  point  in  the  £  coordinate. 
Following  the  traditional  Green’s  function  theory,  we  seek  a  solution  such  that: 


aG[e]  +  bG[4]  +  cG[2]  +  dG  =  S(C  ~  x),  (C.9) 


with  boundary  conditions  such  that  all  boundary  terms  in  (C.8)  are  eliminated: 

G[n]  ( ±il)n  G,  <^[nl  ->  (±il)n  h  n  =  1, 5,  x  ±00,  (C.10) 


where  £(£  —  a:)  is  the  Dirac  delta  function  at  £  =  x.  We  are  interested  in  solutions 
where  a^O.  By  integrating  (C.9)  just  across  £  =  x,  from  x  —  e  to  x  +  e  (e  — ►  0), 
and  requiring  continuity  of  G,G^,G^,G^,  G^,  we  are  left  with: 

aG^\X+t  =  1  (C.ll) 

lx— e 
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Away  from  (  =  iwe  can  write: 


G[6]  +  GiGW  +  a2G[2]  +  a3G  =  0,  (C.12) 

where  oi  =  b/a,  a2  =  c/a ,  03  =  d/a.  Seeking  a  solution  of  the  form: 

G  ~  e,v*  (C.13) 

we  obtain  the  characteristic  polynomial: 

cr6  —  ai<74  —  a2<72  —  a3  =  0.  (C.14) 

For  the  case  in  which  we  are  interested,  the  roots  of  (C.13)  can  be  written  as: 

<7!  =  =  (C.15) 

<72  =  —<75  =  iLi,  (C.16) 

<73  =  <76  =  iL2.  (C.17) 


where  /,  Lx,  L2  are  positive  real  numbers,  and  can  be  obtained  from  the  roots  of 
the  bi-cubic  polynomial  (C.14).  We  now  write  the  solution  for  the  source  function: 

f  G+  =  AGjl{i~x)  +  BGeL *«"*)  +  Coc^W"*)  if  £  <  x 

c?(e,*)  =  ^  +  (c.i8) 

(  G_  =  AGeil(*-V  +  BGeL'(*-V  +  CGeL>(x~V  if  £  >  x 

Continuity  of  G ,  G^,  G^  are  satisfied  automatically,  as  are  the  boundary  condi¬ 
tions  at  ±00.  Continuity  of  G^,  G®,  and  substitution  of  (C.18)  into  (C.ll)  gives 
3  equations  for  the  3  unknowns  AG ,  BG,  and  Cg,  the  solution  being: 


ag  = 

— z 

(C.19) 

2 al  [P  +  Li2)  (/2  +  L22) 

bg  = 

1 

(C.20) 

2 al  (/2  +  Lx2)  [Li2  -  L22) 

CG  = 

1 

(0.21) 

2 al  (P  +  L22)  ( L2 2  -  Li2) 
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(C.19)  can  be  rearranged  to  give: 

Ag  =  2{2 Te 

From  (C.8)  we  can  write: 

Hx)  =  f  G(t,x)gf(t)dt 

J  —  oo 

/X  A  /*-j- oo  A 

G-(i,x)gf(i)dt+  G+((,x)gmd(.  (C.23) 

-OO  J  X 

We  arbitrarily  choose: 

f(x)  =  Ds  exp  (-/?sx2)  .  (C.24) 

For  sufficiently  large  values  of  x  (progressive  wave  traveling  to  greater  values  of 
a:),  and  using  (C.22): 

k*)  =  r  G-((,x)gf(()dt 

J  — oo 

=  gDs  [Aaheilx  +  BGI2e~L'x  +  CGhe~L'x]  ,  (C.25) 

(C.26) 

(C.27) 

(C.28) 
(C.29) 

J2  and  I3  become  negligibly  small  as  x  —>  00,  so: 

j>{x)  «  gDsAGheilx.  (C.30) 


where 


J  exp  [—(3sx2  —  ilx'j  dx  =  ^ 
J  exp  (—/3sx2  +  Lix')  dx  = 
J  exp  {^—fdsx2  +  L2x}  dx  = 


r 

f— exp 


.JL) 

4  (3a)' 

m 

4  M' 


rxp  w. 
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We  here  are  interested  in  waves  propagating  in  the  x  direction.  The  desired 
progressive  wave  solution  (waves  propagating  in  x)  of  (C.l)  and  (C.2)  away  from 
the  source  region  ( x  — >  oo)  is: 


3 

1 

& 

II 

(C.31) 

3 

1 

II 

(C.32) 

7  mo 

9°  w[l  +  C3(khy  +  C,(kh)4} 

(C.33) 

2  1  +  C,(«.)2  +  C2(kh)> 

9  1  +  Cz(kh)2  +  C^kh)4  ’ 

(C.34) 

Setting  A  =  0  (no  y  dependence)  and  l  =  k  we  can  write: 

j>(x,y,t)  =  gDsAGI1ei(kx-“tK  (C.35) 


Substitution  of  (C.32)  and  (C.33)  into  (C.35),  gives  the  relationship  between  the 
source  function  amplitude  Ds  and  the  desired  wave  amplitude  ry0: 


Ds 


_ *?7o _ 

uAGl^l  +  C^khf  +  C^khYY 


(C.36) 
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